Introduction to R (Intensive)

Chris Krogslund
Political Science // D-Lab // UC Berkeley

Where we're going...

Where we're going...

R and Object-oriented Programming (OOP)

The R and RStudio Environments

Getting Help

Arithmetic and Boolean Operations

4 + 5
4 - 5
4 * 5
4/5

log(4)
exp(4)
sin(4)
cos(4)
tan(4)
sqrt(4)

max(4, 5, 6)
min(4, 5, 6)
sum(4, 5, 6)
prod(4, 5, 6)
4==5
4!=5
!4==5

4>5
4!>5 # Caution!
!4>5

4>=5
4!>=5 # Caution
!4>=5

4<5
4!<5 # Caution
!4<5

4<=5
4!<=5 # Caution
!4<=5

4 %in% c(4,5,6)
!4 %in% c(4,5,6)
4 %in% c(5,6)
!4 %in% c(5,6)

Assignment Operators

  1. <-
  2. ->
  3. =
mysum <- 4 + 5
simple.division <- 4/5
15myobject<-4+5
myobject<-4+5
$myobject<-4+5
myobject<-4+5
myobject <- 4 + 5
print(MYOBJECT)
print(myobject)

Assignment Operators/Concatenation

mynumbers <- c(4, 5, 6)
print(mynumbers)
## [1] 4 5 6
4:10
## [1]  4  5  6  7  8  9 10
4:100
##  [1]   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20
## [18]  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37
## [35]  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54
## [52]  55  56  57  58  59  60  61  62  63  64  65  66  67  68  69  70  71
## [69]  72  73  74  75  76  77  78  79  80  81  82  83  84  85  86  87  88
## [86]  89  90  91  92  93  94  95  96  97  98  99 100
seq(0, 100, 10)
##  [1]   0  10  20  30  40  50  60  70  80  90 100
seq(0, 100, 5)
##  [1]   0   5  10  15  20  25  30  35  40  45  50  55  60  65  70  75  80
## [18]  85  90  95 100
rep(6, times = 10)
##  [1] 6 6 6 6 6 6 6 6 6 6
rep(c(1:5), times = 10)
##  [1] 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5
## [36] 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5
rep(c(1:5), each = 10)
##  [1] 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 4 4 4 4 4
## [36] 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5

Object Classes

  1. Numeric
  2. Character
  3. Factor
  4. Logical
  5. Missing/Non-numeric/Null
x <- 4 + 5
class(x)
is.numeric(x)

y <- "New York"
class(y)
is.character(y)
is.numeric(y)

z <- factor(x = c("New York", "New Jersey"), ordered = T)
print(z)
class(z)
is.factor(z)
is.character(z)

a <- TRUE
is.logical(a)
is.character(a)
is.factor(a)

b <- NA
class(b)
is.na(b)
is.numeric(b)

c <- 0/0
class(c)
print(c)
is.nan(c)
is.na(c)

d <- 5/0
class(d)
is.infinite(d)

e <- NULL
class(e)
is.null(e)
is.numeric(e)

Object Classes/Data Structures

  1. Atomic Vectors -- a single piece of data
x <- 3
print(x)
y <- "New York"
print(y)
  1. Vectors -- a uni-dimensional ordered collection of class-homogenous data
x <- c(3, 4, 9, 10)
print(x)
y <- c("New York", "New Jersey", "New Mexico")
print(y)
z <- c(1, 2, "a", "b")
print(z)
class(z)
as.numeric(z)
  1. Matrices -- a two-dimensional ordered collection of class-homogenous data
x <- matrix(data = 1:100, nrow = 10)
print(x)
y <- matrix(data = c(1, 2, "a", "b"), nrow = 2)
print(y)
  1. Arrays -- a generalized n-dimensional ordered collection of class-homogenous data
x <- array(data = 1:300, dim = c(10, 10, 3))
print(x)
  1. Data Frames -- a two-dimensional ordered collection of class-heterogenous data
x <- data.frame(c(1, 2), c("a", "b"))
print(x)
x <- data.frame(col1 = c(1, 2), col2 = c("a", "b"))
print(x)
class(x)
class(x[, 1])
class(x[, 2])
  1. Lists -- a generalized n-dimensional ordered collection of class-heterogenous data
x <- list(3, c(3, 4, "a", "b"), matrix(data = 1:100, nrow = 10), array(data = 1:300, 
    dim = c(10, 10, 3)), data.frame(col1 = c(1, 2), col2 = c("a", "b")))
print(x)
x <- list(myatomic.vector = 3, the.vector = c(3, 4, "a", "b"), some.matrix = matrix(data = 1:100, 
    nrow = 10), another.array = array(data = 1:300, dim = c(10, 10, 3)), onemore.data.frame = data.frame(col1 = c(1, 
    2), col2 = c("a", "b")))
print(x)
class(x)
class(x[[2]])
class(x[[3]])
class(x[[4]])
class(x[[5]])

Programming/Functions/Architecture

  1. Name
  2. Arguments
  3. Content/Body
  4. Return
function.name <- function(argument1, argument2) {
    
    "Content/Body"
    "e.g."
    "object<-argument1*(1-argument2)"
    
    return(object)
    
}

Example #1: Missing Letters

missing.letters1 <- function(text) {
    
    data.split <- strsplit(x = text, split = "")  # Split the text into letters
    data.split <- unlist(data.split)  # Convert the above result from a list to a vector
    data.split <- tolower(data.split)  # Convert all letter to lowercase
    data.split <- unique(data.split)  # Keep only the unique letters
    data.split <- data.split[data.split != " "]  # Get rid of spaces
    missing <- letters[!letters %in% data.split]  # Find the letters that are NOT in the text 
    
    return(missing)  # Return the missing letters
}

missing.letters1(text = "All we have to fear is fear itself")
##  [1] "b" "c" "d" "g" "j" "k" "m" "n" "p" "q" "u" "x" "y" "z"
missing.letters1(text = "Supercalifragilisticexpialidocious")
##  [1] "b" "h" "j" "k" "m" "n" "q" "v" "w" "y" "z"

Or, slightly more compressed:

missing.letters2 <- function(text) {
    
    data.split <- unique(tolower(unlist(strsplit(x = text, split = ""))))
    missing <- letters[!letters %in% data.split[data.split != " "]]
    
    return(missing)
}

identical(missing.letters1(text = "some text"), missing.letters2(text = "some text"))
## [1] TRUE

Programming/Functions/Multi-argument Return

missing.letters2 <- function(text) {
    
    data.split <- unique(tolower(unlist(strsplit(x = text, split = ""))))
    missing <- letters[!letters %in% data.split[data.split != " "]]
    
    num.unique.letters <- length(data.split)
    
    return(missing, num.unique.letters)
}

missing.letters2(text = "Supercalifragilisticexpialidocious")
## Error: multi-argument returns are not permitted

For example:

missing.letters2 <- function(text) {
    
    data.split <- unique(tolower(unlist(strsplit(x = text, split = ""))))
    missing <- letters[!letters %in% data.split[data.split != " "]]
    
    num.unique.letters <- length(data.split)
    
    mylist <- list(missing, num.unique.letters)
    
    return(mylist)
}

missing.letters2(text = "Supercalifragilisticexpialidocious")
## [[1]]
##  [1] "b" "h" "j" "k" "m" "n" "q" "v" "w" "y" "z"
## 
## [[2]]
## [1] 15
example.list <- missing.letters2(text = "Supercalifragilisticexpialidocious")
# Get the first element of the list
example.list[[1]]
##  [1] "b" "h" "j" "k" "m" "n" "q" "v" "w" "y" "z"
# Get the second element of the list
example.list[[2]]
## [1] 15

Programming/Functions/Conditional Execution

if (condition) true.action else false.action
x <- 10

if (x < 10) "fine" else "also fine"
## [1] "also fine"
if (x >= 10) x * 100 else 50/x
## [1] 1000
if (x < 10) {
    y <- 2 * x
    c("really great", y)
} else {
    y <- -2 * x
    c("really not great", y)
}
## [1] "really not great" "-20"
x <- 8:12

if (x < 10) "fine" else "also fine"
## Warning: the condition has length > 1 and only the first element will be
## used
## [1] "fine"
if (x >= 10) x * 100 else 50/x
## Warning: the condition has length > 1 and only the first element will be
## used
## [1] 6.250 5.556 5.000 4.545 4.167
if (x < 10) {
    y <- 2 * x
    c("really great", y)
} else {
    y <- -2 * x
    c("really not great", y)
}
## Warning: the condition has length > 1 and only the first element will be
## used
## [1] "really great" "16"           "18"           "20"          
## [5] "22"           "24"

Programming/Functions/Conditional Execution/ifelse()

Basic usage:

ifelse(test = , yes = , no = )
x = 8:12

ifelse(test = x < 10, yes = "fine", no = "also fine")
## [1] "fine"      "fine"      "also fine" "also fine" "also fine"
ifelse(test = x >= 10, yes = x * 100, no = 50/x)
## [1]    6.250    5.556 1000.000 1100.000 1200.000
ifelse(test = x < 10, yes = {
    y <- 2 * x
    3 * y
}, no = {
    y <- -2 * x
    3 * y
})
## [1]  48  54 -60 -66 -72

Programming/Functions/Conditional Execution/for loops

Basic usage:

storage.object <- NULL
for (element in object.containing.elements) {
    storage.object[element.index] <- operation(element)
}

For example [note a) the vectorization, and b) the storage]

out <- NULL
for (i in 1:length(x)) {
    if (x[i] < 10) {
        y <- 2 * x[i]
        out[i] <- 3 * y
    } else {
        y <- -2 * x[i]
        out[i] <- 3 * y
    }
}
print(out)
## [1]  48  54 -60 -66 -72

Programming/Functions/Conditional Execution/VERSUS

x <- 1:1e+05

time.ifelse <- system.time(ifelse(test = x < 10, yes = {
    y <- 2 * x
    3 * y
}, no = {
    y <- -2 * x
    3 * y
}))
out <- NULL
time.forloop <- system.time(for (i in 1:length(x)) {
    cat(i, "\n")
    if (x[i] < 10) {
        y <- 2 * x[i]
        out[i] <- 3 * y
    } else {
        y <- -2 * x[i]
        out[i] <- 3 * y
    }
})

time.forloop/time.ifelse

Useful Stuff: Applied Research Edition

  1. Importing different types of data from different sources
  2. Summarizing the structure and content of data
  3. Carrying out operations and calculations across groups
  4. Reshaping data to and from various formats
  5. Attempting to conduct causal inference

Importing data/spreadsheets

# Import data with header row, values separated by ',', decimals as '.'
data <- read.csv(file = "  ", stringsAsFactors = )  # Basic Call
# Example
cpds.data <- read.csv(file = "datasets/cpds_comma.csv", stringsAsFactors = F)
dim(cpds.data)  # Check dimensionality of dataset
## [1] 1173  237
cpds.data[1:5, 1:5]  # View first 5 rows/columns
##   year   country countryn gov_right gov_cent
## 1 1960 Australia        1       100        0
## 2 1961 Australia        1       100        0
## 3 1962 Australia        1       100        0
## 4 1963 Australia        1       100        0
## 5 1964 Australia        1       100        0
# Import data with header row, values separated by ';', decimals as ','
data <- read.csv2(file = "  ", stringsAsFactors = )

# Import data with header row, values separated by tab, decimals as '.'
data <- read.delim(file = "  ", stringsAsFactors = )

# Import data with header row, values separated by tab, decimals as ';'
data <- read.delim2(file = "  ", stringsAsFactors = )

# For importing tabular data with maximum customizeability
data <- read.table(file = , header = , sep = , quote = , dec = , fill = , stringsAsFactors = )
# Install the 'gdata' package (only necessary one time)
# install.packages('gdata') # Not Run

# Load the 'gdata' package (necessary every new R session)
library(gdata)
# For importing both .xls and .xlsx files
data <- read.xls(xls = , sheet = )  # Basic Call
# Example with .xls (single sheet)
cpds.data <- read.xls(xls = "datasets/cpds_excel_old.xls")
cpds.data[1:5, 1:5]
##   year   country countryn gov_right gov_cent
## 1 1960 Australia        1       100        0
## 2 1961 Australia        1       100        0
## 3 1962 Australia        1       100        0
## 4 1963 Australia        1       100        0
## 5 1964 Australia        1       100        0
# Example with .xlsx (single sheet)
cpds.data <- read.xls(xls = "datasets/cpds_excel_new.xlsx")
cpds.data[1:5, 1:5]
##   year   country countryn gov_right gov_cent
## 1 1960 Australia        1       100        0
## 2 1961 Australia        1       100        0
## 3 1962 Australia        1       100        0
## 4 1963 Australia        1       100        0
## 5 1964 Australia        1       100        0

Importing data/proprietaries (e.g.: .dta, .spss, .ssd)

# Install the 'foreign' package (only necessary one time)
# install.packages('foreign') # Not Run

# Load the 'foreign' package (necessary every new R session)
library(foreign)

read.xxxx(file = )  # Basic call, where .xxxx is the imported file's extension/client

# For importing Stata files
data <- read.dta(file = )
# Example
cpds.data <- read.dta(file = "datasets/cpds_stata.dta")
cpds.data[1:5, 1:5]
##   year   country countryn gov_right gov_cent
## 1 1960 Australia        1       100        0
## 2 1961 Australia        1       100        0
## 3 1962 Australia        1       100        0
## 4 1963 Australia        1       100        0
## 5 1964 Australia        1       100        0
# For importing SPSS files
data <- read.spss(file = )

# For importing SAS files
data <- read.ssd(file = )
data <- read.xport(file = )

# For importing Fortran files
data <- read.fortran(file = )

# For importing DBF files
data <- read.dbf(file = )

Importing data/urls (e.g.: http, https, ftp)

# Import crime dataset from Quandl
link <- "http://www.quandl.com/api/v1/datasets/DISASTERCENTER/RATES.csv"
crime.data <- read.csv(file = link)
crime.data[1:25, 1:5]
##          Date Total Violent Property Murder
## 1  2012-12-31  3246   386.9     2859    4.7
## 2  2011-12-31  3292   387.1     2905    4.7
## 3  2010-12-31  3350   404.5     2946    4.8
## 4  2009-12-31  3466   431.9     3036    5.0
## 5  2008-12-31  3669   457.5     3212    5.4
## 6  2007-12-31  3730   466.9     3264    5.6
## 7  2006-12-31  3808   473.6     3334    5.7
## 8  2005-12-31  3900   469.0     3432    5.6
## 9  2004-12-31  3977   463.2     3514    5.5
## 10 2003-12-31  4067   475.8     3591    5.7
## 11 2002-12-31  4125   494.4     3631    5.6
## 12 2001-12-31  4163   504.5     3658    5.6
## 13 2000-12-31  4125   506.5     3618    5.5
## 14 1999-12-31  4266   523.0     3744    5.7
## 15 1998-12-31  4616   566.4     4049    6.3
## 16 1997-12-31  4927   611.0     4316    6.8
## 17 1996-12-31  5088   636.6     4451    7.4
## 18 1995-12-31  5275   684.5     4591    8.2
## 19 1994-12-31  5374   713.6     4660    9.0
## 20 1993-12-31  5484   746.8     4738    9.5
## 21 1992-12-31  5660   757.5     4903    9.3
## 22 1991-12-31  5898   758.1     5140    9.8
## 23 1990-12-31  5820   731.8     5088    9.4
## 24 1989-12-31  5741   663.1     5078    8.7
## 25 1988-12-31  5664   637.2     5027    8.4

Importing data/other (e.g.: html tables, html, xml, json)

# Install the 'XML' package (only necessary one time)
# install.packages('XML') # Not Run

# Load the 'XML' package (necessary every new R session)
library(XML)
data <- readHTMLTable(doc = , header = , which = , stringsAsFactors = )  # Basic Call
link <- "http://clerk.house.gov/evs/2014/ROLL_000.asp"
votes.2014 <- readHTMLTable(doc = link, header = T, which = 1, stringsAsFactors = F)

dim(votes.2014)  # Check data dimensionality
## [1] 99  6
colnames(votes.2014)  # Get column names
## [1] "Roll"              "Date"              "Issue"            
## [4] "Question"          "Result"            "Title/Description"
votes.2014[1:25, 1:5]
##    Roll   Date     Issue
## 1    99  6-Mar H RES 501
## 2    98  5-Mar  H R 2126
## 3    97  5-Mar  H R 4118
## 4    96  5-Mar  H R 4118
## 5    95  5-Mar   H R 938
## 6    94  5-Mar H RES 497
## 7    93  5-Mar H RES 497
## 8    92  4-Mar H RES 488
## 9    91  4-Mar  H R 3370
## 10   90 28-Feb   H R 899
## 11   89 28-Feb   H R 899
## 12   88 28-Feb   H R 899
## 13   87 28-Feb   H R 899
## 14   86 28-Feb   H R 899
## 15   85 27-Feb  H R 3193
## 16   84 27-Feb  H R 3193
## 17   83 27-Feb  H R 3193
## 18   82 27-Feb  H R 3193
## 19   81 27-Feb  H R 3193
## 20   80 27-Feb H RES 492
## 21   79 27-Feb H RES 492
## 22   78 27-Feb  H R 2804
## 23   77 27-Feb  H R 2804
## 24   76 27-Feb  H R 2804
## 25   75 27-Feb  H R 2804
##                                                Question Result
## 1                     On Ordering the Previous Question      P
## 2   On Motion to Suspend the Rules and Pass, as Amended      P
## 3                                            On Passage      P
## 4               On Motion to Recommit with Instructions      F
## 5   On Motion to Suspend the Rules and Pass, as Amended      P
## 6                         On Agreeing to the Resolution      P
## 7                     On Ordering the Previous Question      P
## 8  On Motion to Suspend the Rules and Agree, as Amended      P
## 9   On Motion to Suspend the Rules and Pass, as Amended      P
## 10                                           On Passage      P
## 11              On Motion to Recommit with Instructions      F
## 12                         On Agreeing to the Amendment      F
## 13                         On Agreeing to the Amendment      F
## 14                         On Agreeing to the Amendment      F
## 15                                           On Passage      P
## 16              On Motion to Recommit with Instructions      F
## 17                         On Agreeing to the Amendment      F
## 18                         On Agreeing to the Amendment      A
## 19                         On Agreeing to the Amendment      A
## 20                        On Agreeing to the Resolution      P
## 21                    On Ordering the Previous Question      P
## 22                                           On Passage      P
## 23              On Motion to Recommit with Instructions      F
## 24                         On Agreeing to the Amendment      F
## 25                         On Agreeing to the Amendment      F
# Install the 'RCurl' package (only necessary one time)
# install.packages('RCurl') # Not Run

# Load the 'RCurl' and 'XML' packages (necessary every new R session)
library(RCurl)
library(XML)

link <- "http://www.whitehouse.gov/the_press_office/President_Barack_Obamas_Inaugural_Address"
speech <- getURL(url = link)

# Optional html parsing (clean file, extract content)
speech.parsed <- htmlParse(file = speech)

text <- xpathSApply(doc = speech.parsed, path = "//div[@id='content']", fun = xmlValue)

# Optional content analysis

# Install the 'qdap' package (only necessary one time)
# install.packages('qdap') # Not Run

# Load the 'qdap' package (necessary every new R session)
library(qdap)

# Convert the text to a data.frame
text.df <- data.frame(text = text, stringsAsFactors = F)

# Split the text into sentences
sentences <- sentSplit(dataframe = text.df, text.var = "text")

# Compute various polarity (sentiment) statistics
polarity.score <- polarity(text.var = sentences$text)
polarity.score$group[, c("total.sentences", "total.words", "ave.polarity")]
##   total.sentences total.words ave.polarity
## 1             119        2574      0.08423
polarity.score$all[1:20, c("pos.words", "neg.words")]
##                  pos.words                          neg.words
## 1                        -                                  -
## 2                        -                                  -
## 3                        -                                  -
## 4                     free                                  -
## 5          grateful, trust                                  -
## 6  thank, well, generosity                                  -
## 7                        -                                  -
## 8        prosperity, peace                                  -
## 9                        -                             raging
## 10         skill, faithful                                  -
## 11                       -                                  -
## 12                    well                             crisis
## 13                       -                             hatred
## 14                       -        badly, greed, failure, hard
## 15                       -                               lost
## 16                       -             costly, fail, threaten
## 17                       -                             crisis
## 18    profound, confidence nagging, fear, decline, inevitable
## 19                       -                                  -
## 20                       -                                  -
# Specify a link and issue a GET request
link <- "http://rss.nytimes.com/services/xml/rss/nyt/HomePage.xml"
page <- getURL(url = link)

# Parse the results of the get request as an XML file
rss.feed <- xmlParse(file = page)  # Not printed
# Install the 'jsonlite' package (only necessary one time)
# install.packages('jsonlite') # Not Run

# Load the 'jsonlite' package (necessary every new R session)
library(jsonlite)

# Import JSON file, convert to list
cpds.data <- fromJSON(txt = "datasets/cpds_json.txt", simplifyDataFrame = F)
head(cpds.data)
## [[1]]
## [[1]]$year
## [1] 1960
## 
## [[1]]$country
## [1] "Australia"
## 
## [[1]]$variable
## [1] "gov_right"
## 
## [[1]]$value
## [1] 100
## 
## 
## [[2]]
## [[2]]$year
## [1] 1961
## 
## [[2]]$country
## [1] "Australia"
## 
## [[2]]$variable
## [1] "gov_right"
## 
## [[2]]$value
## [1] 100
## 
## 
## [[3]]
## [[3]]$year
## [1] 1962
## 
## [[3]]$country
## [1] "Australia"
## 
## [[3]]$variable
## [1] "gov_right"
## 
## [[3]]$value
## [1] 100
## 
## 
## [[4]]
## [[4]]$year
## [1] 1963
## 
## [[4]]$country
## [1] "Australia"
## 
## [[4]]$variable
## [1] "gov_right"
## 
## [[4]]$value
## [1] 100
## 
## 
## [[5]]
## [[5]]$year
## [1] 1964
## 
## [[5]]$country
## [1] "Australia"
## 
## [[5]]$variable
## [1] "gov_right"
## 
## [[5]]$value
## [1] 100
## 
## 
## [[6]]
## [[6]]$year
## [1] 1965
## 
## [[6]]$country
## [1] "Australia"
## 
## [[6]]$variable
## [1] "gov_right"
## 
## [[6]]$value
## [1] 100
# Import JSON file, convert to data.frame
cpds.data <- fromJSON(txt = "datasets/cpds_json.txt", simplifyDataFrame = T)
head(cpds.data)
##   year   country  variable value
## 1 1960 Australia gov_right   100
## 2 1961 Australia gov_right   100
## 3 1962 Australia gov_right   100
## 4 1963 Australia gov_right   100
## 5 1964 Australia gov_right   100
## 6 1965 Australia gov_right   100

Summarizing data/example dataset

Replication data from Gelman et al., "Rich State, Poor State, Red State, Blue State: What’s the Matter with Connecticut?", Quarterly Journal of Political Science, 2007, 2: 345-367

For decades, the Democrats have been viewed as the party of the poor, with the Republicans representing the rich. Recent presidential elections, however, have shown a reverse pattern, with Democrats performing well in the richer blue states in the northeast and coasts, and Republicans dominating in the red states in the middle of the country and the south. Through multilevel modeling of individual- level survey data and county- and state-level demographic and electoral data, we reconcile these patterns.Furthermore, we find that income matters more in red America than in blue America. In poor states, rich people are much more likely than poor people to vote for the Republican presidential candidate, but in rich states (such as Connecticut), income has a very low correlation with vote preference.

library(foreign)
red.blue <- read.dta("datasets/2004_labeled_processed_race.dta")

Summarizing data/object structure

# Get the object class
class(red.blue)
## [1] "data.frame"
# Get the object dimensionality
dim(red.blue)  # Note this is rows by columns
## [1] 76205    17
# Get the column names
colnames(red.blue)
##  [1] "state"    "pres04"   "sex"      "race"     "age9"     "partyid" 
##  [7] "income"   "relign8"  "age60"    "age65"    "geocode"  "sizeplac"
## [13] "brnagain" "attend"   "year"     "region"   "y"
# Get the row names
rownames(red.blue)[1:50]  # Only the first 50 rows
##  [1] "1"  "2"  "3"  "4"  "5"  "6"  "7"  "8"  "9"  "10" "11" "12" "13" "14"
## [15] "15" "16" "17" "18" "19" "20" "21" "22" "23" "24" "25" "26" "27" "28"
## [29] "29" "30" "31" "32" "33" "34" "35" "36" "37" "38" "39" "40" "41" "42"
## [43] "43" "44" "45" "46" "47" "48" "49" "50"
# View first six rows and all columns
head(red.blue)
##   state pres04    sex  race  age9 partyid income relign8 age60 age65
## 1     2      1 female white 25-29    <NA>   <NA>    <NA> 18-29 25-29
## 2     2      2   male white 18-24    <NA>   <NA>    <NA> 18-29 18-24
## 3     2      1 female black 30-39    <NA>   <NA>    <NA> 30-44 30-39
## 4     2      1 female black 30-39    <NA>   <NA>    <NA> 30-44 30-39
## 5     2      1 female white 40-44    <NA>   <NA>    <NA> 30-44 40-49
## 6     2      1 female white 30-39    <NA>   <NA>    <NA> 30-44 30-39
##   geocode sizeplac brnagain attend year region y
## 1       3    rural     <NA>   <NA> 2004      4 0
## 2       3    rural     <NA>   <NA> 2004      4 1
## 3       3    rural     <NA>   <NA> 2004      4 0
## 4       3    rural     <NA>   <NA> 2004      4 0
## 5       3    rural     <NA>   <NA> 2004      4 0
## 6       3    rural     <NA>   <NA> 2004      4 0
# View last six rows and all columns
tail(red.blue)
##       state pres04    sex  race  age9        partyid          income
## 76200    50      2 female white 65-74     republican   under $15,000
## 76201    50      2 female white 50-59    independent   under $15,000
## 76202    50      2   male white 18-24    independent   under $15,000
## 76203    50      3 female black 30-39 something else   under $15,000
## 76204    50      2 female white 50-59       democrat $30,000-$49,999
## 76205    50      2   male white 45-49     republican $15,000-$29,999
##               relign8      age60      age65 geocode               sizeplac
## 76200        catholic 60 or over 65 or over       1 city: 10,000 to 49,999
## 76201      protestant      45-59      50-64       1                  rural
## 76202 other christian      18-29      18-24       1 city: 10,000 to 49,999
## 76203 other christian      30-44      30-39       1 city: 10,000 to 49,999
## 76204      mormon/lds      45-59      50-64       1 city: 10,000 to 49,999
## 76205  something else      45-59      40-49       2                  rural
##       brnagain attend year region  y
## 76200     <NA>   <NA> 2004      4  1
## 76201     <NA>   <NA> 2004      4  1
## 76202     <NA>   <NA> 2004      4  1
## 76203     <NA>   <NA> 2004      4 NA
## 76204     <NA>   <NA> 2004      4  1
## 76205     <NA>   <NA> 2004      4  1
# Get detailed column-by-column information
str(red.blue)
## 'data.frame':    76205 obs. of  17 variables:
##  $ state   : int  2 2 2 2 2 2 2 2 2 2 ...
##  $ pres04  : int  1 2 1 1 1 1 1 2 2 2 ...
##  $ sex     : Factor w/ 2 levels "male","female": 2 1 2 2 2 2 1 2 2 2 ...
##  $ race    : Factor w/ 5 levels "white","black",..: 1 1 2 2 1 1 1 1 1 1 ...
##  $ age9    : Factor w/ 9 levels "18-24","25-29",..: 2 1 3 3 4 3 4 1 2 1 ...
##  $ partyid : Factor w/ 4 levels "democrat","republican",..: NA NA NA NA NA NA NA NA NA NA ...
##  $ income  : Factor w/ 8 levels "under $15,000",..: NA NA NA NA NA NA NA NA NA NA ...
##  $ relign8 : Factor w/ 8 levels "protestant","catholic",..: NA NA NA NA NA NA NA NA NA NA ...
##  $ age60   : Factor w/ 4 levels "18-29","30-44",..: 1 1 2 2 2 2 2 1 1 1 ...
##  $ age65   : Factor w/ 6 levels "18-24","25-29",..: 2 1 3 3 4 3 4 1 2 1 ...
##  $ geocode : int  3 3 3 3 3 3 3 3 3 3 ...
##  $ sizeplac: Factor w/ 5 levels "city over 500,000",..: 5 5 5 5 5 5 5 5 5 5 ...
##  $ brnagain: Factor w/ 2 levels "yes","no": NA NA NA NA NA NA NA NA NA NA ...
##  $ attend  : Factor w/ 6 levels "more than once a week",..: NA NA NA NA NA NA NA NA NA NA ...
##  $ year    : num  2004 2004 2004 2004 2004 ...
##  $ region  : num  4 4 4 4 4 4 4 4 4 4 ...
##  $ y       : num  0 1 0 0 0 0 0 1 1 1 ...
##  - attr(*, "datalabel")= chr ""
##  - attr(*, "time.stamp")= chr " 6 Jun 2007 14:53"
##  - attr(*, "formats")= chr  "%8.0g" "%8.0g" "%8.0g" "%8.0g" ...
##  - attr(*, "types")= int  251 251 251 251 251 251 251 251 251 251 ...
##  - attr(*, "val.labels")= chr  "stanum" "presak04" "sex" "race" ...
##  - attr(*, "var.labels")= chr  "state id" "in today's election for president, did you just vote for:" "are you:" "are you:" ...
##  - attr(*, "version")= int 8
##  - attr(*, "label.table")=List of 14
##   ..$ stanum  : Named int 2
##   .. ..- attr(*, "names")= chr "alaska"
##   ..$ presak04: Named int  0 1 2 3 9
##   .. ..- attr(*, "names")= chr  "did not vote" "kerry" "bush" "nader" ...
##   ..$ sex     : Named int  1 2
##   .. ..- attr(*, "names")= chr  "male" "female"
##   ..$ race    : Named int  1 2 3 4 5
##   .. ..- attr(*, "names")= chr  "white" "black" "hispanic/latino" "asian" ...
##   ..$ age9    : Named int  1 2 3 4 5 6 7 8 9
##   .. ..- attr(*, "names")= chr  "18-24" "25-29" "30-39" "40-44" ...
##   ..$ partyid : Named int  1 2 3 4
##   .. ..- attr(*, "names")= chr  "democrat" "republican" "independent" "something else"
##   ..$ income  : Named int  1 2 3 4 5 6 7 8
##   .. ..- attr(*, "names")= chr  "under $15,000" "$15,000-$29,999" "$30,000-$49,999" "$50,000-$74,999" ...
##   ..$ relign8 : Named int  1 2 3 4 5 6 7 8
##   .. ..- attr(*, "names")= chr  "protestant" "catholic" "mormon/lds" "other christian" ...
##   ..$ age60   : Named int  1 2 3 4
##   .. ..- attr(*, "names")= chr  "18-29" "30-44" "45-59" "60 or over"
##   ..$ age65   : Named int  1 2 3 4 5 6
##   .. ..- attr(*, "names")= chr  "18-24" "25-29" "30-39" "40-49" ...
##   ..$ geocode : Named int  1 2 3
##   .. ..- attr(*, "names")= chr  "juneau/fairbanks/rural" "anchorage" "anchorage-fairbanks corridor"
##   ..$ sizeplac: Named int  1 2 3 4 5
##   .. ..- attr(*, "names")= chr  "city over 500,000" "city: 50,000 to 500,000" "suburbs" "city: 10,000 to 49,999" ...
##   ..$ brnagain: Named int  1 2
##   .. ..- attr(*, "names")= chr  "yes" "no"
##   ..$ attend  : Named int  1 2 3 4 5 9
##   .. ..- attr(*, "names")= chr  "more than once a week" "once a week" "a few times a month" "a few times a year" ...

Summarizing data/object content

# View unique values observed in 'age9' (note the factor levels)
unique(x = red.blue$age9)
##  [1] 25-29      18-24      30-39      40-44      45-49      65-74     
##  [7] 75 or over 50-59      <NA>       60-64     
## Levels: 18-24 25-29 30-39 40-44 45-49 50-59 60-64 65-74 75 or over
Binned Age True Age
1 18-24
2 25-29
3 30-39
4 40-44
5 45-49
6 50-59
7 60-64
8 65-74
9 75+
  1. Note that these functions are sensitive to missing values (NA); you should be sure to specify na.rm=T to avoid errors
# Sample 100 times from the standard normal distribution
sample.data <- rnorm(n = 100, mean = 0, sd = 1)

# Add some missing values to the sample
sample.data[c(1, 4, 16, 64)] <- NA

# Attempt to calculate the sample mean (presence of NAs)
mean(x = sample.data)
## [1] NA
# Remove missing values from the sample
sample.data <- sample.data[!is.na(sample.data)]

# Attempt to calculate the sample mean (absence of NAs)
mean(x = sample.data)
## [1] -0.1291
  1. These functions are also sensitive to the presence of factor variables; remove the factor levels to avoid errors (usually use one of as.vector(), as.character(), or as.numeric())
# Check for factor levels in 'age9'
is.factor(x = red.blue$age9)
## [1] TRUE
# Attempt to calculate the sample mean ('age9' is factor)
mean(x = red.blue$age9, na.rm = T)
## Warning: argument is not numeric or logical: returning NA
## [1] NA
# Remove factor levels in 'age9'
red.blue$age9 <- as.numeric(x = red.blue$age9)

# Check that there are no more factor levels in 'age9'
is.factor(x = red.blue$age9)
## [1] FALSE
# Attempt to calculate the sample mean ('age9' is not factor)
mean(x = red.blue$age9, na.rm = T)
## [1] 4.544
# Mean
mean(x = as.numeric(red.blue$age9), na.rm = T)
## [1] 4.544
# Median
median(x = as.numeric(red.blue$age9), na.rm = T)
## [1] 4
# Standard Deviation
sd(x = as.numeric(red.blue$age9), na.rm = T)
## [1] 2.238
# Quartiles
quantile(x = as.numeric(red.blue$age9), na.rm = T, probs = seq(from = 0, to = 1, 
    by = 0.25))
##   0%  25%  50%  75% 100% 
##    1    3    4    6    9
# Quintiles
quantile(x = as.numeric(red.blue$age9), na.rm = T, probs = seq(from = 0, to = 1, 
    by = 0.2))
##   0%  20%  40%  60%  80% 100% 
##    1    3    4    5    6    9
# Deciles
quantile(x = as.numeric(red.blue$age9), na.rm = T, probs = seq(from = 0, to = 1, 
    by = 0.1))
##   0%  10%  20%  30%  40%  50%  60%  70%  80%  90% 100% 
##    1    1    3    3    4    4    5    6    6    8    9
# Percentiles
quantile(x = as.numeric(red.blue$age9), na.rm = T, probs = seq(from = 0, to = 1, 
    by = 0.01))
##   0%   1%   2%   3%   4%   5%   6%   7%   8%   9%  10%  11%  12%  13%  14% 
##    1    1    1    1    1    1    1    1    1    1    1    2    2    2    2 
##  15%  16%  17%  18%  19%  20%  21%  22%  23%  24%  25%  26%  27%  28%  29% 
##    2    2    2    2    2    3    3    3    3    3    3    3    3    3    3 
##  30%  31%  32%  33%  34%  35%  36%  37%  38%  39%  40%  41%  42%  43%  44% 
##    3    3    3    3    3    3    3    3    3    4    4    4    4    4    4 
##  45%  46%  47%  48%  49%  50%  51%  52%  53%  54%  55%  56%  57%  58%  59% 
##    4    4    4    4    4    4    5    5    5    5    5    5    5    5    5 
##  60%  61%  62%  63%  64%  65%  66%  67%  68%  69%  70%  71%  72%  73%  74% 
##    5    5    6    6    6    6    6    6    6    6    6    6    6    6    6 
##  75%  76%  77%  78%  79%  80%  81%  82%  83%  84%  85%  86%  87%  88%  89% 
##    6    6    6    6    6    6    6    7    7    7    7    7    7    8    8 
##  90%  91%  92%  93%  94%  95%  96%  97%  98%  99% 100% 
##    8    8    8    8    8    8    9    9    9    9    9
# Compute standard summary statistics for object 'red.blue'
summary(object = red.blue)
##      state          pres04        sex                     race      
##  Min.   : 1.0   Min.   :0.0   male  :34236   white          :60567  
##  1st Qu.:13.0   1st Qu.:1.0   female:41632   black          : 7166  
##  Median :25.0   Median :2.0   NA's  :  337   hispanic/latino: 4430  
##  Mean   :25.7   Mean   :1.6                  asian          : 1028  
##  3rd Qu.:38.0   3rd Qu.:2.0                  other          : 1781  
##  Max.   :50.0   Max.   :9.0                  NA's           : 1233  
##                 NA's   :543                                         
##       age9               partyid                  income     
##  Min.   :1.0   democrat      :27059   $50,000-$74,999:18864  
##  1st Qu.:3.0   republican    :25940   under $15,000  :16416  
##  Median :4.0   independent   :16512   $30,000-$49,999:15745  
##  Mean   :4.5   something else: 3409   $15,000-$29,999:15563  
##  3rd Qu.:6.0   NA's          : 3285   $75,000-$99,999: 2274  
##  Max.   :9.0                          (Other)        :    0  
##  NA's   :505                          NA's           : 7343  
##             relign8             age60              age65      
##  protestant     :21783   18-29     :14704   18-24     : 8256  
##  catholic       :16411   30-44     :23159   25-29     : 6448  
##  other christian:15230   45-59     :23742   30-39     :14212  
##  none           : 7540   60 or over:14095   40-49     :17938  
##  something else : 4412   NA's      :  505   50-64     :19746  
##  (Other)        : 3794                      65 or over: 9100  
##  NA's           : 7035                      NA's      :  505  
##     geocode                        sizeplac     brnagain    
##  Min.   :1.00   city over 500,000      : 5882   yes :14034  
##  1st Qu.:1.00   city: 50,000 to 500,000:15462   no  :16647  
##  Median :2.00   suburbs                :28796   NA's:45524  
##  Mean   :2.52   city: 10,000 to 49,999 : 8449               
##  3rd Qu.:3.00   rural                  :17501               
##  Max.   :6.00   NA's                   :  115               
##                                                             
##                    attend           year          region    
##  more than once a week: 3669   Min.   :2004   Min.   :1.00  
##  once a week          : 7255   1st Qu.:2004   1st Qu.:1.00  
##  a few times a month  : 4146   Median :2004   Median :3.00  
##  a few times a year   : 8232   Mean   :2004   Mean   :2.45  
##  never                : 4549   3rd Qu.:2004   3rd Qu.:4.00  
##  omit                 :  687   Max.   :2004   Max.   :4.00  
##  NA's                 :47667                                
##        y       
##  Min.   :0.0   
##  1st Qu.:0.0   
##  Median :0.0   
##  Mean   :0.5   
##  3rd Qu.:1.0   
##  Max.   :1.0   
##  NA's   :1787
# Isolate the gender column, remove factor levels
gender <- as.character(x = red.blue$sex)

# Remove missing values
gender <- gender[!is.na(gender)]

# View unique values
unique(gender)
## [1] "female" "male"
# Subset the gender vector to include only males
males <- gender[gender == "male"]
head(males)
## [1] "male" "male" "male" "male" "male" "male"
# Compute the percentage of males using vector lengths
length(males)
## [1] 34236
length(gender)
## [1] 75868
length(males)/length(gender) * 100
## [1] 45.13

Group-wise Operations

All techniques for this problem rely on the split-apply-combine strategy

First, take the data (or some object) and split it into smaller datasets on the basis of some variable

Dataset A

x y z
1 1 1
2 2 1
3 3 1
4 1 2
5 2 2
6 3 2

Datasets B and C (Dataset A split according to "z")

x y z x y z
1 1 1 4 1 2
2 2 1 5 2 2
3 3 1 6 3 2

Second, apply some function to each one of the smaller datasets/objects

Example function: mean of variables "x" and "y"

Datasets B' and C'

mean(x) mean(y) z mean(x) mean(y) z
2 2 1 5 2 2

Third, combine the results into a larger dataset/object

Datasets B' and C'

mean(x) mean(y) z mean(x) mean(y) z
2 2 1 5 2 2

Dataset A'

mean(x) mean(y) z
2 2 1
5 2 2

Group-wise Operations/plyr

  1. Don't have to worry about different name, argument, or output consistencies
  2. Easily parallelized
  3. Input from, and output to, data frames, matricies, and lists
  4. Progress bars for lengthy computation
  5. Informative error messages
# Install the 'plyr' package (only necessary one time)
# install.packages('plyr') # Not Run

# Load the 'plyr' package (necessary every new R session)
library(plyr)

Group-wise Operations/plyr/selecting functions

  1. What is the class of your input object?
  2. What is the class of your desired output object?

Group-wise Operations/plyr/writing commands

All of the major plyr functions have the same basic syntax

xxply(.data = , .variables = , .fun = , ...)

Consider the case where we want to calculate vote choice statistics across race from a data.frame, and return them as a data.frame:

# Check for (and remove) factor levels in 'race' and 'pres04'
is.factor(x = red.blue$race)
## [1] TRUE
red.blue$race <- as.character(x = red.blue$race)
is.factor(x = red.blue$race)
## [1] FALSE
is.factor(x = red.blue$pres04)
## [1] FALSE
# Check the class of the object
class(red.blue)
## [1] "data.frame"
# Using the appropriate plyr function (ddply), compute vote percentages for
# Kerry (pres04==1), Bush (pres04==2), Nader (pres04==3), and others
# (pres04==9)
ddply(.data = red.blue, .variables = .(race), .fun = summarize, kerry = length(pres04[pres04 == 
    1])/length(pres04) * 100, bush = length(pres04[pres04 == 2])/length(pres04) * 
    100, nader = length(pres04[pres04 == 3])/length(pres04) * 100, other = length(pres04[pres04 == 
    9])/length(pres04) * 100)
##              race kerry  bush nader  other
## 1           asian 61.67 38.13 1.459 0.9728
## 2           black 87.18 12.39 1.675 1.1862
## 3 hispanic/latino 60.93 37.77 1.535 1.8736
## 4           other 59.40 37.90 2.471 3.0320
## 5           white 43.87 55.20 1.327 1.3142
## 6            <NA> 55.07 45.09 2.920 3.0008

Consider the case where we want to calculate vote choice statistics across race from a data.frame, and return them as a list:

dlply(.data = red.blue, .variables = .(race), .fun = summarize, kerry = length(pres04[pres04 == 
    1])/length(pres04) * 100, bush = length(pres04[pres04 == 2])/length(pres04) * 
    100, nader = length(pres04[pres04 == 3])/length(pres04) * 100, other = length(pres04[pres04 == 
    9])/length(pres04) * 100)
## $asian
##   kerry  bush nader  other
## 1 61.67 38.13 1.459 0.9728
## 
## $black
##   kerry  bush nader other
## 1 87.18 12.39 1.675 1.186
## 
## $`hispanic/latino`
##   kerry  bush nader other
## 1 60.93 37.77 1.535 1.874
## 
## $other
##   kerry bush nader other
## 1  59.4 37.9 2.471 3.032
## 
## $white
##   kerry bush nader other
## 1 43.87 55.2 1.327 1.314
## 
## $`NA`
##   kerry  bush nader other
## 1 55.07 45.09  2.92 3.001
## 
## attr(,"split_type")
## [1] "data.frame"
## attr(,"split_labels")
##              race
## 1           asian
## 2           black
## 3 hispanic/latino
## 4           other
## 5           white
## 6            <NA>

Consider the case where we want to calculate vote choice statistics across race from a list, and return them as a data.frame:

# Split the data.frame into a list on the basis of 'race'
red.blue.race.list <- split(x = red.blue, f = red.blue$race)

# Check the class of the object
class(red.blue.race.list)
## [1] "list"
# Check for list element names
objects(red.blue.race.list)
## [1] "asian"           "black"           "hispanic/latino" "other"          
## [5] "white"
# Compute summary statistics (note: no .variable argument)
ldply(.data = red.blue.race.list, .fun = summarize, kerry = length(pres04[pres04 == 
    1])/length(pres04) * 100, bush = length(pres04[pres04 == 2])/length(pres04) * 
    100, nader = length(pres04[pres04 == 3])/length(pres04) * 100, other = length(pres04[pres04 == 
    9])/length(pres04) * 100)
##               .id kerry  bush nader  other
## 1           asian 61.67 38.13 1.459 0.9728
## 2           black 87.18 12.39 1.675 1.1862
## 3 hispanic/latino 60.93 37.77 1.535 1.8736
## 4           other 59.40 37.90 2.471 3.0320
## 5           white 43.87 55.20 1.327 1.3142

Consider the case where we want to calculate vote choice statistics across race from a list, and return them as a list:

llply(.data = red.blue.race.list, .fun = summarize, kerry = length(pres04[pres04 == 
    1])/length(pres04) * 100, bush = length(pres04[pres04 == 2])/length(pres04) * 
    100, nader = length(pres04[pres04 == 3])/length(pres04) * 100, other = length(pres04[pres04 == 
    9])/length(pres04) * 100)
## $asian
##   kerry  bush nader  other
## 1 61.67 38.13 1.459 0.9728
## 
## $black
##   kerry  bush nader other
## 1 87.18 12.39 1.675 1.186
## 
## $`hispanic/latino`
##   kerry  bush nader other
## 1 60.93 37.77 1.535 1.874
## 
## $other
##   kerry bush nader other
## 1  59.4 37.9 2.471 3.032
## 
## $white
##   kerry bush nader other
## 1 43.87 55.2 1.327 1.314

Group-wise Operations/plyr/functions

transform: applies a function to a data.frame and adds new vectors (columns) to it

# Add a column containing the average age of the race of the individual
red.blue.transform <- ddply(.data = red.blue, .variables = .(race), .fun = transform, 
    race.avg.age = mean(x = age9, na.rm = T))
unique(red.blue.transform$race.avg.age)
## [1] 4.133 4.003 3.704 4.029 4.674 5.464

Note that transform can't do transformations that involve the results of other transformations from the same call

# Attempt to add new columns that draw on other (but still new) columns
red.blue.transform <- ddply(.data = red.blue, .variables = .(race), .fun = transform, 
    race.avg.age = mean(x = age9, na.rm = T), race.avg.age.plusone = race.avg.age + 
        1)
## Error: object 'race.avg.age' not found

For this, we need mutate: just like transform, but it executes the commands iteratively so transformations can be carried out that rely on previous transformations from the same call

# Attempt to add new columns that draw on other (but still new) columns
red.blue.mutate <- ddply(.data = red.blue, .variables = .(race), .fun = mutate, 
    race.avg.age = mean(x = age9, na.rm = T), race.avg.age.plusone = race.avg.age + 
        1)
unique(red.blue.mutate$race.avg.age)
## [1] 4.133 4.003 3.704 4.029 4.674 5.464
unique(red.blue.mutate$race.avg.age.plusone)
## [1] 5.133 5.003 4.704 5.029 5.674 6.464

Another very useful function is arrange, which orders a data frame on the basis of column contents

# Arrange by 'age9', ascending
red.blue.age <- plyr::arrange(df = red.blue, age9)
red.blue.age[1:25, 1:5]
##    state pres04    sex            race age9
## 1      2      2   male           white    1
## 2      2      2 female           white    1
## 3      2      2 female           white    1
## 4      2      2   male           black    1
## 5      2      2   male           white    1
## 6      2      1   male           white    1
## 7      2      1 female hispanic/latino    1
## 8      2      1   male           black    1
## 9      2      1   male           white    1
## 10     2      1   male           white    1
## 11     2      1 female           white    1
## 12     2      2   male hispanic/latino    1
## 13     2      2   male           white    1
## 14     2      1   male           white    1
## 15     2      1 female           white    1
## 16     2      1 female           white    1
## 17     2      2 female           white    1
## 18     2      2   male           other    1
## 19     2      1 female hispanic/latino    1
## 20     2     NA female hispanic/latino    1
## 21     2      2   male           white    1
## 22     2      1 female           white    1
## 23     2      2 female           white    1
## 24     2      1 female           white    1
## 25     2      1   male           white    1
# Arrange by 'age9', descending
red.blue.age <- plyr::arrange(df = red.blue, desc(age9))
red.blue.age[1:25, 1:5]
##    state pres04    sex  race age9
## 1      2      9   male white    9
## 2      2      1 female white    9
## 3      2      2   male white    9
## 4      2      1   male white    9
## 5      2      2 female white    9
## 6      2      1   male white    9
## 7      2      2 female white    9
## 8      2      1 female white    9
## 9      2      2 female white    9
## 10     2      2   male white    9
## 11     2      2 female white    9
## 12     2      1 female white    9
## 13     2      2 female asian    9
## 14     2      2   male white    9
## 15     2      2   male white    9
## 16     2      1 female white    9
## 17     1      2   male white    9
## 18     1      1   male white    9
## 19     1      1 female  <NA>    9
## 20     1      2   male white    9
## 21     1      2 female white    9
## 22     1      2   male white    9
## 23     1      1   male white    9
## 24     1      2 female white    9
## 25     1     NA female white    9
# Arrange by 'age9' then 'sex'
red.blue.age.sex <- plyr::arrange(df = red.blue, age9, sex)
red.blue.age.sex[1:25, 1:5]
##    state pres04  sex            race age9
## 1      2      2 male           white    1
## 2      2      2 male           black    1
## 3      2      2 male           white    1
## 4      2      1 male           white    1
## 5      2      1 male           black    1
## 6      2      1 male           white    1
## 7      2      1 male           white    1
## 8      2      2 male hispanic/latino    1
## 9      2      2 male           white    1
## 10     2      1 male           white    1
## 11     2      2 male           other    1
## 12     2      2 male           white    1
## 13     2      1 male           white    1
## 14     2      2 male           white    1
## 15     2      2 male           white    1
## 16     2      2 male           white    1
## 17     2      2 male hispanic/latino    1
## 18     2      1 male           white    1
## 19     2      2 male           white    1
## 20     2      1 male           white    1
## 21     2      2 male           white    1
## 22     2      2 male           white    1
## 23     2      2 male           white    1
## 24     2      9 male           white    1
## 25     2      2 male           white    1
# Arrange by 'sex' (descending) then 'age9'
red.blue.sex.age <- plyr::arrange(df = red.blue, desc(sex), age9)
red.blue.sex.age[1:25, 1:5]
##    state pres04    sex            race age9
## 1      2      2 female           white    1
## 2      2      2 female           white    1
## 3      2      1 female hispanic/latino    1
## 4      2      1 female           white    1
## 5      2      1 female           white    1
## 6      2      1 female           white    1
## 7      2      2 female           white    1
## 8      2      1 female hispanic/latino    1
## 9      2     NA female hispanic/latino    1
## 10     2      1 female           white    1
## 11     2      2 female           white    1
## 12     2      1 female           white    1
## 13     2      2 female           white    1
## 14     2      1 female           white    1
## 15     2      2 female           white    1
## 16     2      2 female           white    1
## 17     2     NA female           white    1
## 18     2      1 female           white    1
## 19     2      2 female           other    1
## 20     2      1 female           white    1
## 21     2      2 female           white    1
## 22     2      2 female            <NA>    1
## 23     2      2 female           white    1
## 24     2      3 female           white    1
## 25     2      9 female           white    1

Group-wise Operations/dplyr

# CAUTION: this is a really big file
load(file = "datasets/ushouse_votes")

# View first six rows, check dimensionality
head(ushouse.votes)
##           name party state    vote year vote.num yea.vote
## 1     Ackerman     D    NY Present 1990        1        0
## 2        Akaka     D    HI Present 1990        1        0
## 3    Alexander     D    AR Present 1990        1        0
## 4     Anderson     D    CA Present 1990        1        0
## 5 Andrews (TX)     D    TX Present 1990        1        0
## 6     Annunzio     D    IL Present 1990        1        0
dim(ushouse.votes)
## [1] 6754366       7
# Install the 'dplyr' package (only necessary one time)
# install.packages('dplyr') # Not Run

# Load the 'dplyr' package (necessary every new R session)
library(dplyr)
# Basic syntax for a dplyr call:

# 1) Create a tbl (tabular data structure) that groups the data
data <- group_by(x = , ...)

# 2) Execute the pre-defined function
summarise(.data = , ...)

A side-by-side comparison of ddply and dplyr (with benchmarks):

# ddply
time.ddply <- system.time(avg.yea.states.ddply <- ddply(.data = ushouse.votes, 
    .variables = .(state), summarize, avg.yea = mean(x = yea.vote, na.rm = T)))
avg.yea.states.ddply
##    state avg.yea
## 1     AK  0.5346
## 2     AL  0.5960
## 3     AR  0.6135
## 4     AZ  0.6048
## 5     CA  0.6000
## 6     CO  0.6084
## 7     CT  0.6322
## 8     DE  0.6375
## 9     FL  0.6040
## 10    GA  0.6121
## 11    HI  0.6091
## 12    IA  0.6213
## 13    ID  0.6057
## 14    IL  0.6009
## 15    IN  0.6117
## 16    KS  0.6254
## 17    KY  0.6160
## 18    LA  0.5890
## 19    MA  0.6127
## 20    MD  0.6274
## 21    ME  0.6342
## 22    MI  0.6111
## 23    MN  0.6135
## 24    MO  0.5985
## 25    MS  0.6116
## 26    MT  0.5923
## 27    NC  0.6222
## 28    ND  0.6412
## 29    NE  0.6233
## 30    NH  0.6544
## 31    NJ  0.6156
## 32    NM  0.6170
## 33    NV  0.6329
## 34    NY  0.6098
## 35    OH  0.6139
## 36    OK  0.6017
## 37    OR  0.6122
## 38    PA  0.6092
## 39    RI  0.6160
## 40    SC  0.6193
## 41    SD  0.6641
## 42    TN  0.6229
## 43    TX  0.5958
## 44    UT  0.6163
## 45    VA  0.6125
## 46    VT  0.6072
## 47    WA  0.6161
## 48    WI  0.6307
## 49    WV  0.6040
## 50    WY  0.5404
## 51    XX  0.3112
# dplyr
state.groups <- group_by(x = ushouse.votes, state)
time.dplyr <- system.time(avg.yea.states.dplyr <- summarise(.data = state.groups, 
    avg.yea = mean(x = yea.vote, na.rm = T)))

# Compare the execution times (ratio of ddply-to-plyr)
time.ddply/time.dplyr
##    user  system elapsed 
##   16.50   18.01   18.31

Reshaping Data/reshape2

# Install the 'reshape2' package (only necessary one time)
# install.packages('reshape2') # Not Run

# Load the 'reshape2' package (necessary every new R session)
library(reshape2)

Reshaping Data/reshape2/melt

  1. Unit-of-Analysis identifiers, or columns you don't want to melt down
  2. Characteristic variables, or columns you do want to melt down
# Basic Call
melt(data = , id.vars = , measure.vars = , variable.name = , value.name = )

To see how this works in practice, consider a subset of the "red.blue" data.frame containing only the first 10 rows and 5 columns

red.blue.subset <- red.blue[1:10, 1:5]
red.blue.subset
##    state pres04    sex  race age9
## 1      2      1 female white    2
## 2      2      2   male white    1
## 3      2      1 female black    3
## 4      2      1 female black    3
## 5      2      1 female white    4
## 6      2      1 female white    3
## 7      2      1   male white    4
## 8      2      2 female white    1
## 9      2      2 female white    2
## 10     2      2 female white    1

Suppose we wanted to convert this data from its current wide format to an entirely long format. How to proceed?

First, select which columns you want to keep (i.e. not melt). In this case, I'm interested in having individual voters as my unit of analysis. Unfortunately, there is no column containing an individual identification number in this data, so I'll just add one as "individual":

red.blue.subset$individual <- 1:nrow(red.blue.subset)
red.blue.subset
##    state pres04    sex  race age9 individual
## 1      2      1 female white    2          1
## 2      2      2   male white    1          2
## 3      2      1 female black    3          3
## 4      2      1 female black    3          4
## 5      2      1 female white    4          5
## 6      2      1 female white    3          6
## 7      2      1   male white    4          7
## 8      2      2 female white    1          8
## 9      2      2 female white    2          9
## 10     2      2 female white    1         10

Second, select which columns you want to melt. In this case, I'd like to melt every column except "individual".

With these two pieces of information, I'm ready to melt down the data.frame:

melt(data = red.blue.subset, id.vars = "individual", measure.vars = c("state", 
    "pres04", "sex", "race", "age9"))
##    individual variable  value
## 1           1    state      2
## 2           2    state      2
## 3           3    state      2
## 4           4    state      2
## 5           5    state      2
## 6           6    state      2
## 7           7    state      2
## 8           8    state      2
## 9           9    state      2
## 10         10    state      2
## 11          1   pres04      1
## 12          2   pres04      2
## 13          3   pres04      1
## 14          4   pres04      1
## 15          5   pres04      1
## 16          6   pres04      1
## 17          7   pres04      1
## 18          8   pres04      2
## 19          9   pres04      2
## 20         10   pres04      2
## 21          1      sex female
## 22          2      sex   male
## 23          3      sex female
## 24          4      sex female
## 25          5      sex female
## 26          6      sex female
## 27          7      sex   male
## 28          8      sex female
## 29          9      sex female
## 30         10      sex female
## 31          1     race  white
## 32          2     race  white
## 33          3     race  black
## 34          4     race  black
## 35          5     race  white
## 36          6     race  white
## 37          7     race  white
## 38          8     race  white
## 39          9     race  white
## 40         10     race  white
## 41          1     age9      2
## 42          2     age9      1
## 43          3     age9      3
## 44          4     age9      3
## 45          5     age9      4
## 46          6     age9      3
## 47          7     age9      4
## 48          8     age9      1
## 49          9     age9      2
## 50         10     age9      1
# If you want to melt ALL columns that aren't ID variables, you can also
# omit the 'measure.vars' argument
melt(data = red.blue.subset, id.vars = "individual")
##    individual variable  value
## 1           1    state      2
## 2           2    state      2
## 3           3    state      2
## 4           4    state      2
## 5           5    state      2
## 6           6    state      2
## 7           7    state      2
## 8           8    state      2
## 9           9    state      2
## 10         10    state      2
## 11          1   pres04      1
## 12          2   pres04      2
## 13          3   pres04      1
## 14          4   pres04      1
## 15          5   pres04      1
## 16          6   pres04      1
## 17          7   pres04      1
## 18          8   pres04      2
## 19          9   pres04      2
## 20         10   pres04      2
## 21          1      sex female
## 22          2      sex   male
## 23          3      sex female
## 24          4      sex female
## 25          5      sex female
## 26          6      sex female
## 27          7      sex   male
## 28          8      sex female
## 29          9      sex female
## 30         10      sex female
## 31          1     race  white
## 32          2     race  white
## 33          3     race  black
## 34          4     race  black
## 35          5     race  white
## 36          6     race  white
## 37          7     race  white
## 38          8     race  white
## 39          9     race  white
## 40         10     race  white
## 41          1     age9      2
## 42          2     age9      1
## 43          3     age9      3
## 44          4     age9      3
## 45          5     age9      4
## 46          6     age9      3
## 47          7     age9      4
## 48          8     age9      1
## 49          9     age9      2
## 50         10     age9      1

Note that melt collapses all of the measure variables into two columns: one containing the column/measurement name, the other containing the column/measurement value for that row. By default, these columns are named "variable" and "value", though they can be customized using the "variable.name" and "value.name" arguments. For example:

melt(data = red.blue.subset, id.vars = "individual", measure.vars = c("state", 
    "pres04", "sex", "race", "age9"), variable.name = "characteristic", value.name = "response")
##    individual characteristic response
## 1           1          state        2
## 2           2          state        2
## 3           3          state        2
## 4           4          state        2
## 5           5          state        2
## 6           6          state        2
## 7           7          state        2
## 8           8          state        2
## 9           9          state        2
## 10         10          state        2
## 11          1         pres04        1
## 12          2         pres04        2
## 13          3         pres04        1
## 14          4         pres04        1
## 15          5         pres04        1
## 16          6         pres04        1
## 17          7         pres04        1
## 18          8         pres04        2
## 19          9         pres04        2
## 20         10         pres04        2
## 21          1            sex   female
## 22          2            sex     male
## 23          3            sex   female
## 24          4            sex   female
## 25          5            sex   female
## 26          6            sex   female
## 27          7            sex     male
## 28          8            sex   female
## 29          9            sex   female
## 30         10            sex   female
## 31          1           race    white
## 32          2           race    white
## 33          3           race    black
## 34          4           race    black
## 35          5           race    white
## 36          6           race    white
## 37          7           race    white
## 38          8           race    white
## 39          9           race    white
## 40         10           race    white
## 41          1           age9        2
## 42          2           age9        1
## 43          3           age9        3
## 44          4           age9        3
## 45          5           age9        4
## 46          6           age9        3
## 47          7           age9        4
## 48          8           age9        1
## 49          9           age9        2
## 50         10           age9        1

Note also that one need not melt down all columns that aren't serving as ID columns. The melted data.frame will only contain the values of the measure variables you select. For instance:

red.blue.melt <- melt(data = red.blue.subset, id.vars = "individual", measure.vars = c("pres04", 
    "sex"))
red.blue.melt
##    individual variable  value
## 1           1   pres04      1
## 2           2   pres04      2
## 3           3   pres04      1
## 4           4   pres04      1
## 5           5   pres04      1
## 6           6   pres04      1
## 7           7   pres04      1
## 8           8   pres04      2
## 9           9   pres04      2
## 10         10   pres04      2
## 11          1      sex female
## 12          2      sex   male
## 13          3      sex female
## 14          4      sex female
## 15          5      sex female
## 16          6      sex female
## 17          7      sex   male
## 18          8      sex female
## 19          9      sex female
## 20         10      sex female

Reshaping Data/reshape2/cast

dcast(data = , formula = xvar1 + xvar2 ~ yvar1 + yvar2, value.var = , fun.aggregate = )

Some example usages:

# Original data
red.blue.subset
##    state pres04    sex  race age9 individual
## 1      2      1 female white    2          1
## 2      2      2   male white    1          2
## 3      2      1 female black    3          3
## 4      2      1 female black    3          4
## 5      2      1 female white    4          5
## 6      2      1 female white    3          6
## 7      2      1   male white    4          7
## 8      2      2 female white    1          8
## 9      2      2 female white    2          9
## 10     2      2 female white    1         10
# Cast a data.frame containing the individual column and columns containing
# the expansion of 'age9' on the basis of its unique values
dcast(data = red.blue.subset, formula = individual ~ age9, value.var = "age9")
##    individual  1  2  3  4
## 1           1 NA  2 NA NA
## 2           2  1 NA NA NA
## 3           3 NA NA  3 NA
## 4           4 NA NA  3 NA
## 5           5 NA NA NA  4
## 6           6 NA NA  3 NA
## 7           7 NA NA NA  4
## 8           8  1 NA NA NA
## 9           9 NA  2 NA NA
## 10         10  1 NA NA NA
# Previously melted data
red.blue.melt
##    individual variable  value
## 1           1   pres04      1
## 2           2   pres04      2
## 3           3   pres04      1
## 4           4   pres04      1
## 5           5   pres04      1
## 6           6   pres04      1
## 7           7   pres04      1
## 8           8   pres04      2
## 9           9   pres04      2
## 10         10   pres04      2
## 11          1      sex female
## 12          2      sex   male
## 13          3      sex female
## 14          4      sex female
## 15          5      sex female
## 16          6      sex female
## 17          7      sex   male
## 18          8      sex female
## 19          9      sex female
## 20         10      sex female
# Cast a new data.frame from melted data.frame containing the individual
# column and expanding the 'variable' column
dcast(data = red.blue.melt, formula = individual ~ variable, value.var = "value")
##    individual pres04    sex
## 1           1      1 female
## 2           2      2   male
## 3           3      1 female
## 4           4      1 female
## 5           5      1 female
## 6           6      1 female
## 7           7      1   male
## 8           8      2 female
## 9           9      2 female
## 10         10      2 female

Inference

  1. Dealing with missing values
  2. Carrying out Classical Hypothesis Tests
  3. Estimating Regressions
  4. Carryingout Regression Diagnostics

Inference/Missing Values

# Install the 'Amelia' package (only necessary one time)
# install.packages('Amelia') # Not Run

# Load the 'Amelia' package (necessary every new R session)
library(Amelia)

Consider the following panel dataset containing economic and demographic data for several African countries over roughly two decades starting in the 1970s.

data(africa)  # Load the 'africa' data.set (included in Amelia)

# Check dimensionality, colnames, head, tail, and structure
dim(africa)
## [1] 120   7
colnames(africa)
## [1] "year"       "country"    "gdp_pc"     "infl"       "trade"     
## [6] "civlib"     "population"
head(africa)
##   year      country gdp_pc  infl trade civlib population
## 1 1972 Burkina Faso    377 -2.92 29.69 0.5000    5848380
## 2 1973 Burkina Faso    376  7.60 31.31 0.5000    5958700
## 3 1974 Burkina Faso    393  8.72 35.22 0.3333    6075700
## 4 1975 Burkina Faso    416 18.76 40.11 0.3333    6202000
## 5 1976 Burkina Faso    435 -8.40 37.76 0.5000    6341030
## 6 1977 Burkina Faso    448 29.99 41.11 0.6667    6486870
tail(africa)
##     year country gdp_pc   infl trade civlib population
## 115 1986  Zambia    786  51.85 86.02 0.3333    6901040
## 116 1987  Zambia    760  43.01 75.16 0.3333    7107800
## 117 1988  Zambia    757  55.56 59.47 0.3333    7322250
## 118 1989  Zambia    756 127.89 60.63 0.3333    7546740
## 119 1990  Zambia    689 117.51 72.47 0.3333    7784000
## 120 1991  Zambia    699  92.60 71.86 0.6667    8022380
str(africa)
## 'data.frame':    120 obs. of  7 variables:
##  $ year      : num  1972 1973 1974 1975 1976 ...
##  $ country   : Factor w/ 6 levels "Burkina Faso",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ gdp_pc    : int  377 376 393 416 435 448 445 461 457 466 ...
##  $ infl      : num  -2.92 7.6 8.72 18.76 -8.4 ...
##  $ trade     : num  29.7 31.3 35.2 40.1 37.8 ...
##  $ civlib    : num  0.5 0.5 0.333 0.333 0.5 ...
##  $ population: int  5848380 5958700 6075700 6202000 6341030 6486870 6639120 6797540 6962000 7132320 ...

We can check for missing values in each of the variables by calling summary()

summary(africa)  # Note the missing values in 'gdp_pc' and 'trade'
##       year              country       gdp_pc          infl       
##  Min.   :1972   Burkina Faso:20   Min.   : 376   Min.   : -8.40  
##  1st Qu.:1977   Burundi     :20   1st Qu.: 514   1st Qu.:  4.76  
##  Median :1982   Cameroon    :20   Median :1036   Median :  8.72  
##  Mean   :1982   Congo       :20   Mean   :1058   Mean   : 12.75  
##  3rd Qu.:1986   Senegal     :20   3rd Qu.:1245   3rd Qu.: 13.56  
##  Max.   :1991   Zambia      :20   Max.   :2723   Max.   :127.89  
##                                   NA's   :2                      
##      trade           civlib        population      
##  Min.   : 24.4   Min.   :0.000   Min.   : 1332490  
##  1st Qu.: 38.5   1st Qu.:0.167   1st Qu.: 4332190  
##  Median : 59.6   Median :0.167   Median : 5853565  
##  Mean   : 62.6   Mean   :0.289   Mean   : 5765594  
##  3rd Qu.: 81.2   3rd Qu.:0.333   3rd Qu.: 7355000  
##  Max.   :134.1   Max.   :0.667   Max.   :11825390  
##  NA's   :5

A distribution of potential values for these missing values can be imputed using the Amelia package.

# Basic Call
amelia(x = , m = , ts = , cs = )

For the "africa" dataset:

# Run 5 imputations
africa.imputed <- amelia(x = africa, m = 5, ts = "year", cs = "country")
## -- Imputation 1 --
## 
##   1  2  3
## 
## -- Imputation 2 --
## 
##   1  2
## 
## -- Imputation 3 --
## 
##   1  2  3
## 
## -- Imputation 4 --
## 
##   1  2  3
## 
## -- Imputation 5 --
## 
##   1  2
# Note the class of the resulting object
class(africa.imputed)
## [1] "amelia"
# Note the many objects it contains
objects(africa.imputed)
##  [1] "arguments"   "code"        "covMatrices" "imputations" "iterHist"   
##  [6] "m"           "message"     "missMatrix"  "mu"          "orig.vars"  
## [11] "overvalues"  "theta"
# Most important object is likely 'imputations' (this contains a list of
# data.frames with imputed values included).  Must be sure to apply your
# analyses (e.g. regressions) over each data.frame in the list, and then
# aggregate the results
str(africa.imputed$imputations)
## List of 5
##  $ imp1:'data.frame':    120 obs. of  7 variables:
##   ..$ year      : num [1:120] 1972 1973 1974 1975 1976 ...
##   ..$ country   : Factor w/ 6 levels "Burkina Faso",..: 1 1 1 1 1 1 1 1 1 1 ...
##   ..$ gdp_pc    : num [1:120] 377 376 393 416 435 448 445 461 457 466 ...
##   ..$ infl      : num [1:120] -2.92 7.6 8.72 18.76 -8.4 ...
##   ..$ trade     : num [1:120] 29.7 31.3 35.2 40.1 37.8 ...
##   ..$ civlib    : num [1:120] 0.5 0.5 0.333 0.333 0.5 ...
##   ..$ population: int [1:120] 5848380 5958700 6075700 6202000 6341030 6486870 6639120 6797540 6962000 7132320 ...
##  $ imp2:'data.frame':    120 obs. of  7 variables:
##   ..$ year      : num [1:120] 1972 1973 1974 1975 1976 ...
##   ..$ country   : Factor w/ 6 levels "Burkina Faso",..: 1 1 1 1 1 1 1 1 1 1 ...
##   ..$ gdp_pc    : num [1:120] 377 376 393 416 435 448 445 461 457 466 ...
##   ..$ infl      : num [1:120] -2.92 7.6 8.72 18.76 -8.4 ...
##   ..$ trade     : num [1:120] 29.7 31.3 35.2 40.1 37.8 ...
##   ..$ civlib    : num [1:120] 0.5 0.5 0.333 0.333 0.5 ...
##   ..$ population: int [1:120] 5848380 5958700 6075700 6202000 6341030 6486870 6639120 6797540 6962000 7132320 ...
##  $ imp3:'data.frame':    120 obs. of  7 variables:
##   ..$ year      : num [1:120] 1972 1973 1974 1975 1976 ...
##   ..$ country   : Factor w/ 6 levels "Burkina Faso",..: 1 1 1 1 1 1 1 1 1 1 ...
##   ..$ gdp_pc    : num [1:120] 377 376 393 416 435 448 445 461 457 466 ...
##   ..$ infl      : num [1:120] -2.92 7.6 8.72 18.76 -8.4 ...
##   ..$ trade     : num [1:120] 29.7 31.3 35.2 40.1 37.8 ...
##   ..$ civlib    : num [1:120] 0.5 0.5 0.333 0.333 0.5 ...
##   ..$ population: int [1:120] 5848380 5958700 6075700 6202000 6341030 6486870 6639120 6797540 6962000 7132320 ...
##  $ imp4:'data.frame':    120 obs. of  7 variables:
##   ..$ year      : num [1:120] 1972 1973 1974 1975 1976 ...
##   ..$ country   : Factor w/ 6 levels "Burkina Faso",..: 1 1 1 1 1 1 1 1 1 1 ...
##   ..$ gdp_pc    : num [1:120] 377 376 393 416 435 448 445 461 457 466 ...
##   ..$ infl      : num [1:120] -2.92 7.6 8.72 18.76 -8.4 ...
##   ..$ trade     : num [1:120] 29.7 31.3 35.2 40.1 37.8 ...
##   ..$ civlib    : num [1:120] 0.5 0.5 0.333 0.333 0.5 ...
##   ..$ population: int [1:120] 5848380 5958700 6075700 6202000 6341030 6486870 6639120 6797540 6962000 7132320 ...
##  $ imp5:'data.frame':    120 obs. of  7 variables:
##   ..$ year      : num [1:120] 1972 1973 1974 1975 1976 ...
##   ..$ country   : Factor w/ 6 levels "Burkina Faso",..: 1 1 1 1 1 1 1 1 1 1 ...
##   ..$ gdp_pc    : num [1:120] 377 376 393 416 435 448 445 461 457 466 ...
##   ..$ infl      : num [1:120] -2.92 7.6 8.72 18.76 -8.4 ...
##   ..$ trade     : num [1:120] 29.7 31.3 35.2 40.1 37.8 ...
##   ..$ civlib    : num [1:120] 0.5 0.5 0.333 0.333 0.5 ...
##   ..$ population: int [1:120] 5848380 5958700 6075700 6202000 6341030 6486870 6639120 6797540 6962000 7132320 ...
##  - attr(*, "class")= chr [1:2] "mi" "list"
# Can also apply summary methods to objects of class 'amelia'
summary(africa.imputed)
## 
## Amelia output with 5 imputed datasets.
## Return code:  1 
## Message:  Normal EM convergence. 
## 
## Chain Lengths:
## --------------
## Imputation 1:  3
## Imputation 2:  2
## Imputation 3:  3
## Imputation 4:  3
## Imputation 5:  2
## 
## Rows after Listwise Deletion:  115 
## Rows after Imputation:  120 
## Patterns of missingness in the data:  3 
## 
## Fraction Missing for original variables: 
## -----------------------------------------
## 
##            Fraction Missing
## year                0.00000
## country             0.00000
## gdp_pc              0.01667
## infl                0.00000
## trade               0.04167
## civlib              0.00000
## population          0.00000

An example of carring out regressions over each of the imputed datasets:

# Extract the imputed datasets from the amelia object
datasets <- africa.imputed$imputations
# Apply the regression over each element of the list
lapply(X = datasets, FUN = function(x) {
    summary(lm(formula = gdp_pc ~ infl + trade, data = x))$coefficients
})
## $imp1
##             Estimate Std. Error  t value  Pr(>|t|)
## (Intercept)    1.128     80.966  0.01394 9.889e-01
## infl          -6.884      1.656 -4.15596 6.194e-05
## trade         18.466      1.205 15.33040 9.278e-30
## 
## $imp2
##             Estimate Std. Error t value  Pr(>|t|)
## (Intercept)   11.571     80.324  0.1441 8.857e-01
## infl          -6.927      1.648 -4.2026 5.183e-05
## trade         18.377      1.197 15.3518 8.318e-30
## 
## $imp3
##             Estimate Std. Error  t value  Pr(>|t|)
## (Intercept)   -4.544     82.390 -0.05515 9.561e-01
## infl          -6.842      1.659 -4.12359 7.005e-05
## trade         18.507      1.221 15.16127 2.208e-29
## 
## $imp4
##             Estimate Std. Error t value  Pr(>|t|)
## (Intercept)    7.681     83.767  0.0917 9.271e-01
## infl          -6.860      1.681 -4.0816 8.210e-05
## trade         18.374      1.243 14.7805 1.573e-28
## 
## $imp5
##             Estimate Std. Error  t value  Pr(>|t|)
## (Intercept)   -5.878     81.753 -0.07191 9.428e-01
## infl          -6.922      1.656 -4.17887 5.676e-05
## trade         18.584      1.212 15.33062 9.267e-30

Inference/Hypothesis Tests

Suppose we have two different samples, A and B, both drawn from the standard normal distribution:

a <- rnorm(n = 5000, mean = 0, sd = 1)
b <- rnorm(n = 5000, mean = 0, sd = 1)

Suppose we also have a third sample, C, drawn from the normal distribution with mean=1 and sd=0:

c <- rnorm(n = 5000, mean = 1, sd = 1)

One can test for differences in these distributions in either a) their means using t-tests, or b) their entire distributions using ks-tests

# Basic Call
t.test(x = , y = , var.equal = , conf.level = , formula = )
ks.test(x = , y = )
t.test(x = a, y = b)  # No difference in means
## 
##  Welch Two Sample t-test
## 
## data:  a and b
## t = 1.295, df = 9997, p-value = 0.1953
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.01315  0.06434
## sample estimates:
## mean of x mean of y 
##  0.019855 -0.005742
t.test(x = a, y = c)  # Difference in means
## 
##  Welch Two Sample t-test
## 
## data:  a and c
## t = -49.69, df = 9997, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -1.0299 -0.9517
## sample estimates:
## mean of x mean of y 
##   0.01985   1.01067
ks.test(x = a, y = b)  # No difference in distributions
## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  a and b
## D = 0.0186, p-value = 0.3527
## alternative hypothesis: two-sided
ks.test(x = a, y = c)  # Difference in distributions
## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  a and c
## D = 0.3852, p-value < 2.2e-16
## alternative hypothesis: two-sided

Inference/Regression

glm(formula = y ~ x1 + x2 + x3 + ..., family = familyname(link = "linkname"), 
    data = )
# Create an indicator for men aged 65 or over
red.blue$oldman <- ifelse(test = (red.blue$age65 == "65 or over" & red.blue$sex == 
    "male"), yes = 1, no = 0)

# Regress being a 65+ year old male on partyid, income, and religion (logit
# model)
oldman.reg <- glm(formula = oldman ~ partyid + income + relign8, family = binomial(link = "logit"), 
    data = red.blue)
# View objects contained in the regression output
objects(oldman.reg)
##  [1] "aic"               "boundary"          "call"             
##  [4] "coefficients"      "contrasts"         "control"          
##  [7] "converged"         "data"              "deviance"         
## [10] "df.null"           "df.residual"       "effects"          
## [13] "family"            "fitted.values"     "formula"          
## [16] "iter"              "linear.predictors" "method"           
## [19] "model"             "na.action"         "null.deviance"    
## [22] "offset"            "prior.weights"     "qr"               
## [25] "R"                 "rank"              "residuals"        
## [28] "terms"             "weights"           "xlevels"          
## [31] "y"
# Examine regression coefficients
oldman.reg$coefficients
##            (Intercept)      partyidrepublican     partyidindependent 
##                -2.0125                 0.0951                 0.1283 
##  partyidsomething else  income$15,000-$29,999  income$30,000-$49,999 
##                -0.7358                -0.1543                -0.5768 
##  income$50,000-$74,999  income$75,000-$99,999        relign8catholic 
##                -0.8934                -0.8447                -0.5387 
##      relign8mormon/lds relign8other christian          relign8jewish 
##                -0.5416                -1.8544                 0.1297 
##          relign8muslim  relign8something else            relign8none 
##                -1.7141                -1.7263                -1.1830
# Examine regression DoF
oldman.reg$df.residual
## [1] 64191
# Examine regression fit (AIC)
oldman.reg$aic
## [1] 23741
summary(oldman.reg)
## 
## Call:
## glm(formula = oldman ~ partyid + income + relign8, family = binomial(link = "logit"), 
##     data = red.blue)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -0.565  -0.377  -0.269  -0.197   3.220  
## 
## Coefficients:
##                        Estimate Std. Error z value Pr(>|z|)    
## (Intercept)             -2.0125     0.0454  -44.33  < 2e-16 ***
## partyidrepublican        0.0951     0.0441    2.15  0.03118 *  
## partyidindependent       0.1283     0.0492    2.61  0.00914 ** 
## partyidsomething else   -0.7358     0.1361   -5.41  6.4e-08 ***
## income$15,000-$29,999   -0.1543     0.0489   -3.15  0.00161 ** 
## income$30,000-$49,999   -0.5768     0.0533  -10.82  < 2e-16 ***
## income$50,000-$74,999   -0.8934     0.0544  -16.44  < 2e-16 ***
## income$75,000-$99,999   -0.8447     0.1191   -7.09  1.3e-12 ***
## relign8catholic         -0.5387     0.0454  -11.87  < 2e-16 ***
## relign8mormon/lds       -0.5416     0.1051   -5.15  2.5e-07 ***
## relign8other christian  -1.8544     0.0727  -25.50  < 2e-16 ***
## relign8jewish            0.1297     0.1056    1.23  0.21958    
## relign8muslim           -1.7141     0.4524   -3.79  0.00015 ***
## relign8something else   -1.7263     0.1223  -14.11  < 2e-16 ***
## relign8none             -1.1830     0.0775  -15.26  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 25272  on 64205  degrees of freedom
## Residual deviance: 23711  on 64191  degrees of freedom
##   (11999 observations deleted due to missingness)
## AIC: 23741
## 
## Number of Fisher Scoring iterations: 6
# Store summary method results
sum.oldman.reg <- summary(oldman.reg)
# View summary method results objects
objects(sum.oldman.reg)
##  [1] "aic"            "aliased"        "call"           "coefficients"  
##  [5] "contrasts"      "cov.scaled"     "cov.unscaled"   "deviance"      
##  [9] "deviance.resid" "df"             "df.null"        "df.residual"   
## [13] "dispersion"     "family"         "iter"           "na.action"     
## [17] "null.deviance"  "terms"
# View table of coefficients
sum.oldman.reg$coefficients
##                        Estimate Std. Error z value   Pr(>|z|)
## (Intercept)             -2.0125    0.04540 -44.328  0.000e+00
## partyidrepublican        0.0951    0.04414   2.155  3.118e-02
## partyidindependent       0.1283    0.04921   2.607  9.145e-03
## partyidsomething else   -0.7358    0.13607  -5.407  6.394e-08
## income$15,000-$29,999   -0.1543    0.04894  -3.154  1.612e-03
## income$30,000-$49,999   -0.5768    0.05332 -10.817  2.850e-27
## income$50,000-$74,999   -0.8934    0.05436 -16.436  1.062e-60
## income$75,000-$99,999   -0.8447    0.11915  -7.090  1.346e-12
## relign8catholic         -0.5387    0.04540 -11.867  1.765e-32
## relign8mormon/lds       -0.5416    0.10508  -5.154  2.546e-07
## relign8other christian  -1.8544    0.07273 -25.496 2.183e-143
## relign8jewish            0.1297    0.10563   1.228  2.196e-01
## relign8muslim           -1.7141    0.45244  -3.789  1.515e-04
## relign8something else   -1.7263    0.12234 -14.110  3.273e-45
## relign8none             -1.1830    0.07752 -15.261  1.394e-52
glm(formula = y ~ x1 + x2 + factor(x3), family = family(link = "link"), data = )

x1:x2 interacts all terms in x1 with all terms in x2

summary(glm(formula = oldman ~ partyid:income, family = binomial(link = "logit"), 
    data = red.blue))
## 
## Call:
## glm(formula = oldman ~ partyid:income, family = binomial(link = "logit"), 
##     data = red.blue)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -0.398  -0.361  -0.309  -0.271   3.074  
## 
## Coefficients: (1 not defined because of singularities)
##                                             Estimate Std. Error z value
## (Intercept)                                  -2.7213     0.4617   -5.89
## partyiddemocrat:incomeunder $15,000           0.0255     0.4641    0.05
## partyidrepublican:incomeunder $15,000         0.2251     0.4654    0.48
## partyidindependent:incomeunder $15,000        0.0307     0.4669    0.07
## partyidsomething else:incomeunder $15,000    -1.0980     0.5117   -2.15
## partyiddemocrat:income$15,000-$29,999        -0.1651     0.4654   -0.35
## partyidrepublican:income$15,000-$29,999       0.1527     0.4648    0.33
## partyidindependent:income$15,000-$29,999      0.0692     0.4667    0.15
## partyidsomething else:income$15,000-$29,999  -1.0909     0.5263   -2.07
## partyiddemocrat:income$30,000-$49,999        -0.4840     0.4671   -1.04
## partyidrepublican:income$30,000-$49,999      -0.2968     0.4657   -0.64
## partyidindependent:income$30,000-$49,999     -0.2756     0.4683   -0.59
## partyidsomething else:income$30,000-$49,999  -1.0996     0.5349   -2.06
## partyiddemocrat:income$50,000-$74,999        -0.7158     0.4676   -1.53
## partyidrepublican:income$50,000-$74,999      -0.5196     0.4656   -1.12
## partyidindependent:income$50,000-$74,999     -0.5620     0.4687   -1.20
## partyidsomething else:income$50,000-$74,999  -1.9957     0.6175   -3.23
## partyiddemocrat:income$75,000-$99,999        -0.9575     0.5306   -1.80
## partyidrepublican:income$75,000-$99,999      -0.2822     0.4835   -0.58
## partyidindependent:income$75,000-$99,999     -0.6631     0.5270   -1.26
## partyidsomething else:income$75,000-$99,999       NA         NA      NA
##                                             Pr(>|z|)    
## (Intercept)                                  3.8e-09 ***
## partyiddemocrat:incomeunder $15,000           0.9562    
## partyidrepublican:incomeunder $15,000         0.6286    
## partyidindependent:incomeunder $15,000        0.9476    
## partyidsomething else:incomeunder $15,000     0.0319 *  
## partyiddemocrat:income$15,000-$29,999         0.7227    
## partyidrepublican:income$15,000-$29,999       0.7426    
## partyidindependent:income$15,000-$29,999      0.8821    
## partyidsomething else:income$15,000-$29,999   0.0382 *  
## partyiddemocrat:income$30,000-$49,999         0.3001    
## partyidrepublican:income$30,000-$49,999       0.5240    
## partyidindependent:income$30,000-$49,999      0.5562    
## partyidsomething else:income$30,000-$49,999   0.0398 *  
## partyiddemocrat:income$50,000-$74,999         0.1258    
## partyidrepublican:income$50,000-$74,999       0.2644    
## partyidindependent:income$50,000-$74,999      0.2305    
## partyidsomething else:income$50,000-$74,999   0.0012 ** 
## partyiddemocrat:income$75,000-$99,999         0.0711 .  
## partyidrepublican:income$75,000-$99,999       0.5594    
## partyidindependent:income$75,000-$99,999      0.2083    
## partyidsomething else:income$75,000-$99,999       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 26801  on 68088  degrees of freedom
## Residual deviance: 26433  on 68069  degrees of freedom
##   (8116 observations deleted due to missingness)
## AIC: 26473
## 
## Number of Fisher Scoring iterations: 7

**x1/*x2** produces the cross of x1 and x2, or x1+x2+x1:x2

summary(glm(formula = oldman ~ partyid * income, family = binomial(link = "logit"), 
    data = red.blue))
## 
## Call:
## glm(formula = oldman ~ partyid * income, family = binomial(link = "logit"), 
##     data = red.blue)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -0.398  -0.361  -0.309  -0.271   3.074  
## 
## Coefficients:
##                                             Estimate Std. Error z value
## (Intercept)                                 -2.69582    0.04731  -56.99
## partyidrepublican                            0.19963    0.07525    2.65
## partyidindependent                           0.00519    0.08421    0.06
## partyidsomething else                       -1.12346    0.22561   -4.98
## income$15,000-$29,999                       -0.19062    0.07505   -2.54
## income$30,000-$49,999                       -0.50951    0.08509   -5.99
## income$50,000-$74,999                       -0.74125    0.08791   -8.43
## income$75,000-$99,999                       -0.98301    0.26568   -3.70
## partyidrepublican:income$15,000-$29,999      0.11817    0.10922    1.08
## partyidindependent:income$15,000-$29,999     0.22917    0.12313    1.86
## partyidsomething else:income$15,000-$29,999  0.19770    0.34377    0.58
## partyidrepublican:income$30,000-$49,999     -0.01239    0.11999   -0.10
## partyidindependent:income$30,000-$49,999     0.20323    0.13517    1.50
## partyidsomething else:income$30,000-$49,999  0.50794    0.35902    1.41
## partyidrepublican:income$50,000-$74,999     -0.00342    0.12150   -0.03
## partyidindependent:income$50,000-$74,999     0.14855    0.13838    1.07
## partyidsomething else:income$50,000-$74,999 -0.15647    0.47386   -0.33
## partyidrepublican:income$75,000-$99,999      0.47565    0.30756    1.55
## partyidindependent:income$75,000-$99,999     0.28925    0.37425    0.77
## partyidsomething else:income$75,000-$99,999  2.08100    0.57655    3.61
##                                             Pr(>|z|)    
## (Intercept)                                  < 2e-16 ***
## partyidrepublican                            0.00798 ** 
## partyidindependent                           0.95086    
## partyidsomething else                        6.4e-07 ***
## income$15,000-$29,999                        0.01109 *  
## income$30,000-$49,999                        2.1e-09 ***
## income$50,000-$74,999                        < 2e-16 ***
## income$75,000-$99,999                        0.00022 ***
## partyidrepublican:income$15,000-$29,999      0.27930    
## partyidindependent:income$15,000-$29,999     0.06272 .  
## partyidsomething else:income$15,000-$29,999  0.56524    
## partyidrepublican:income$30,000-$49,999      0.91776    
## partyidindependent:income$30,000-$49,999     0.13271    
## partyidsomething else:income$30,000-$49,999  0.15713    
## partyidrepublican:income$50,000-$74,999      0.97756    
## partyidindependent:income$50,000-$74,999     0.28306    
## partyidsomething else:income$50,000-$74,999  0.74124    
## partyidrepublican:income$75,000-$99,999      0.12198    
## partyidindependent:income$75,000-$99,999     0.43959    
## partyidsomething else:income$75,000-$99,999  0.00031 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 26801  on 68088  degrees of freedom
## Residual deviance: 26433  on 68069  degrees of freedom
##   (8116 observations deleted due to missingness)
## AIC: 26473
## 
## Number of Fisher Scoring iterations: 7
# Install the 'devtools' package (only necessary one time)
# install.packages('devtools') # Not Run

# Load the 'devtools' package (necessary every new R session)
library(devtools)

# Source the 'glm.sim' function from gist 9359460
source_gist("https://gist.github.com/ckrogs/9359460")
# Basic Call
glm.sim(data = , dv = , graphs = , graph.best = , robust = , interactions = , 
    quadratics = )

Consider an example with the "africa" dataset. Suppose we wanted to regress "gdp_pc" on some combination of "infl", "trade", "civlib", and "population", but we were unsure which model best explained the outcome. glm.sim can runs all the possible models and aggregates the results.

# Remove extraneous columns from the data
africa.subset <- africa[, c(3:7)]

# Run glm.sim with graphs
output <- glm.sim(data = africa.subset, dv = "gdp_pc", graphs = T)
## Creating model set...
## Number of model terms: 1... 2... 3... 4... Done.
## Running glm...
## Percent complete: 6.67 %
## 13.33 %
## 20 %
## 26.67 %
## 33.33 %
## 40 %
## 46.67 %
## 53.33 %
## 60 %
## 66.67 %
## 73.33 %
## 80 %
## 86.67 %
## 93.33 %
## 100 %
## Done.
## Identifying best models (est)...
## 6.67 %
## 13.33 %
## 20 %
## 26.67 %
## 33.33 %
## 40 %
## 46.67 %
## 53.33 %
## 60 %
## 66.67 %
## 73.33 %
## 80 %
## 86.67 %
## 93.33 %
## 100 %
## Binding data...
## Arranging data...
## Identifying best models (p)...
## 6.67 %
## 13.33 %
## 20 %
## 26.67 %
## 33.33 %
## 40 %
## 46.67 %
## 53.33 %
## 60 %
## 66.67 %
## 73.33 %
## 80 %
## 86.67 %
## 93.33 %
## 100 %
## Binding data...
## Arranging data...
## Done.
## Extracting results...
## 6.67 %
## 13.33 %
## 20 %
## 26.67 %
## 33.33 %
## 40 %
## 46.67 %
## 53.33 %
## 60 %
## 66.67 %
## 73.33 %
## 80 %
## 86.67 %
## 93.33 %
## 100 %
## Done.
## Creating graphs... Done.

# Note the objects
objects(output)
## [1] "best.est"   "best.p"     "coef.table" "plot.est"   "plot.p"

# Run glm.sim with graphs and the best model
output <- glm.sim(data = africa.subset, dv = "gdp_pc", graphs = T, graph.best = 1)
## Creating model set...
## Number of model terms: 1... 2... 3... 4... Done.
## Running glm...
## Percent complete: 6.67 %
## 13.33 %
## 20 %
## 26.67 %
## 33.33 %
## 40 %
## 46.67 %
## 53.33 %
## 60 %
## 66.67 %
## 73.33 %
## 80 %
## 86.67 %
## 93.33 %
## 100 %
## Done.
## Identifying best models (est)...
## 6.67 %
## 13.33 %
## 20 %
## 26.67 %
## 33.33 %
## 40 %
## 46.67 %
## 53.33 %
## 60 %
## 66.67 %
## 73.33 %
## 80 %
## 86.67 %
## 93.33 %
## 100 %
## Binding data...
## Arranging data...
## Identifying best models (p)...
## 6.67 %
## 13.33 %
## 20 %
## 26.67 %
## 33.33 %
## 40 %
## 46.67 %
## 53.33 %
## 60 %
## 66.67 %
## 73.33 %
## 80 %
## 86.67 %
## 93.33 %
## 100 %
## Binding data...
## Arranging data...
## Done.
## Extracting results...
## 6.67 %
## 13.33 %
## 20 %
## 26.67 %
## 33.33 %
## 40 %
## 46.67 %
## 53.33 %
## 60 %
## 66.67 %
## 73.33 %
## 80 %
## 86.67 %
## 93.33 %
## 100 %
## Done.
## Creating graphs... Done.
output$plot.p
output$plot.est
# Run glm.sim with graphs, best model, and robust standard errors
output <- glm.sim(data = africa.subset, dv = "gdp_pc", graphs = T, graph.best = 1, 
    robust = T)
## Creating model set...
## Number of model terms: 1... 2... 3... 4... Done.
## Running glm...
## Percent complete: 6.67 %
## 13.33 %
## 20 %
## 26.67 %
## 33.33 %
## 40 %
## 46.67 %
## 53.33 %
## 60 %
## 66.67 %
## 73.33 %
## 80 %
## 86.67 %
## 93.33 %
## 100 %
## Done.
## Identifying best models (est)...
## 6.67 %
## 13.33 %
## 20 %
## 26.67 %
## 33.33 %
## 40 %
## 46.67 %
## 53.33 %
## 60 %
## 66.67 %
## 73.33 %
## 80 %
## 86.67 %
## 93.33 %
## 100 %
## Binding data...
## Arranging data...
## Identifying best models (p)...
## 6.67 %
## 13.33 %
## 20 %
## 26.67 %
## 33.33 %
## 40 %
## 46.67 %
## 53.33 %
## 60 %
## 66.67 %
## 73.33 %
## 80 %
## 86.67 %
## 93.33 %
## 100 %
## Binding data...
## Arranging data...
## Done.
## Extracting results...
## 6.67 %
## 13.33 %
## 20 %
## 26.67 %
## 33.33 %
## 40 %
## 46.67 %
## 53.33 %
## 60 %
## 66.67 %
## 73.33 %
## 80 %
## 86.67 %
## 93.33 %
## 100 %
## Done.
## Creating graphs... Done.
output$plot.p
output$plot.est
# Run glm.sim with graphs, best model, and first-order interactions
output <- glm.sim(data = africa.subset, dv = "gdp_pc", graphs = T, graph.best = 1, 
    interactions = T)
## Creating model set...
## Number of model terms: 1... 2... 3... 4... 5... 6... 7... 8... 9... 10... Done.
## Running glm...
## Percent complete: 0.1 %
## 0.2 %
## 0.29 %
## 0.39 %
## 0.49 %
## 0.59 %
## 0.68 %
## 0.78 %
## 0.88 %
## 0.98 %
## 1.08 %
## 1.17 %
## 1.27 %
## 1.37 %
## 1.47 %
## 1.56 %
## 1.66 %
## 1.76 %
## 1.86 %
## 1.96 %
## 2.05 %
## 2.15 %
## 2.25 %
## 2.35 %
## 2.44 %
## 2.54 %
## 2.64 %
## 2.74 %
## 2.83 %
## 2.93 %
## 3.03 %
## 3.13 %
## 3.23 %
## 3.32 %
## 3.42 %
## 3.52 %
## 3.62 %
## 3.71 %
## 3.81 %
## 3.91 %
## 4.01 %
## 4.11 %
## 4.2 %
## 4.3 %
## 4.4 %
## 4.5 %
## 4.59 %
## 4.69 %
## 4.79 %
## 4.89 %
## 4.99 %
## 5.08 %
## 5.18 %
## 5.28 %
## 5.38 %
## 5.47 %
## 5.57 %
## 5.67 %
## 5.77 %
## 5.87 %
## 5.96 %
## 6.06 %
## 6.16 %
## 6.26 %
## 6.35 %
## 6.45 %
## 6.55 %
## 6.65 %
## 6.74 %
## 6.84 %
## 6.94 %
## 7.04 %
## 7.14 %
## 7.23 %
## 7.33 %
## 7.43 %
## 7.53 %
## 7.62 %
## 7.72 %
## 7.82 %
## 7.92 %
## 8.02 %
## 8.11 %
## 8.21 %
## 8.31 %
## 8.41 %
## 8.5 %
## 8.6 %
## 8.7 %
## 8.8 %
## 8.9 %
## 8.99 %
## 9.09 %
## 9.19 %
## 9.29 %
## 9.38 %
## 9.48 %
## 9.58 %
## 9.68 %
## 9.78 %
## 9.87 %
## 9.97 %
## 10.07 %
## 10.17 %
## 10.26 %
## 10.36 %
## 10.46 %
## 10.56 %
## 10.65 %
## 10.75 %
## 10.85 %
## 10.95 %
## 11.05 %
## 11.14 %
## 11.24 %
## 11.34 %
## 11.44 %
## 11.53 %
## 11.63 %
## 11.73 %
## 11.83 %
## 11.93 %
## 12.02 %
## 12.12 %
## 12.22 %
## 12.32 %
## 12.41 %
## 12.51 %
## 12.61 %
## 12.71 %
## 12.81 %
## 12.9 %
## 13 %
## 13.1 %
## 13.2 %
## 13.29 %
## 13.39 %
## 13.49 %
## 13.59 %
## 13.69 %
## 13.78 %
## 13.88 %
## 13.98 %
## 14.08 %
## 14.17 %
## 14.27 %
## 14.37 %
## 14.47 %
## 14.57 %
## 14.66 %
## 14.76 %
## 14.86 %
## 14.96 %
## 15.05 %
## 15.15 %
## 15.25 %
## 15.35 %
## 15.44 %
## 15.54 %
## 15.64 %
## 15.74 %
## 15.84 %
## 15.93 %
## 16.03 %
## 16.13 %
## 16.23 %
## 16.32 %
## 16.42 %
## 16.52 %
## 16.62 %
## 16.72 %
## 16.81 %
## 16.91 %
## 17.01 %
## 17.11 %
## 17.2 %
## 17.3 %
## 17.4 %
## 17.5 %
## 17.6 %
## 17.69 %
## 17.79 %
## 17.89 %
## 17.99 %
## 18.08 %
## 18.18 %
## 18.28 %
## 18.38 %
## 18.48 %
## 18.57 %
## 18.67 %
## 18.77 %
## 18.87 %
## 18.96 %
## 19.06 %
## 19.16 %
## 19.26 %
## 19.35 %
## 19.45 %
## 19.55 %
## 19.65 %
## 19.75 %
## 19.84 %
## 19.94 %
## 20.04 %
## 20.14 %
## 20.23 %
## 20.33 %
## 20.43 %
## 20.53 %
## 20.63 %
## 20.72 %
## 20.82 %
## 20.92 %
## 21.02 %
## 21.11 %
## 21.21 %
## 21.31 %
## 21.41 %
## 21.51 %
## 21.6 %
## 21.7 %
## 21.8 %
## 21.9 %
## 21.99 %
## 22.09 %
## 22.19 %
## 22.29 %
## 22.39 %
## 22.48 %
## 22.58 %
## 22.68 %
## 22.78 %
## 22.87 %
## 22.97 %
## 23.07 %
## 23.17 %
## 23.26 %
## 23.36 %
## 23.46 %
## 23.56 %
## 23.66 %
## 23.75 %
## 23.85 %
## 23.95 %
## 24.05 %
## 24.14 %
## 24.24 %
## 24.34 %
## 24.44 %
## 24.54 %
## 24.63 %
## 24.73 %
## 24.83 %
## 24.93 %
## 25.02 %
## 25.12 %
## 25.22 %
## 25.32 %
## 25.42 %
## 25.51 %
## 25.61 %
## 25.71 %
## 25.81 %
## 25.9 %
## 26 %
## 26.1 %
## 26.2 %
## 26.3 %
## 26.39 %
## 26.49 %
## 26.59 %
## 26.69 %
## 26.78 %
## 26.88 %
## 26.98 %
## 27.08 %
## 27.17 %
## 27.27 %
## 27.37 %
## 27.47 %
## 27.57 %
## 27.66 %
## 27.76 %
## 27.86 %
## 27.96 %
## 28.05 %
## 28.15 %
## 28.25 %
## 28.35 %
## 28.45 %
## 28.54 %
## 28.64 %
## 28.74 %
## 28.84 %
## 28.93 %
## 29.03 %
## 29.13 %
## 29.23 %
## 29.33 %
## 29.42 %
## 29.52 %
## 29.62 %
## 29.72 %
## 29.81 %
## 29.91 %
## 30.01 %
## 30.11 %
## 30.21 %
## 30.3 %
## 30.4 %
## 30.5 %
## 30.6 %
## 30.69 %
## 30.79 %
## 30.89 %
## 30.99 %
## 31.09 %
## 31.18 %
## 31.28 %
## 31.38 %
## 31.48 %
## 31.57 %
## 31.67 %
## 31.77 %
## 31.87 %
## 31.96 %
## 32.06 %
## 32.16 %
## 32.26 %
## 32.36 %
## 32.45 %
## 32.55 %
## 32.65 %
## 32.75 %
## 32.84 %
## 32.94 %
## 33.04 %
## 33.14 %
## 33.24 %
## 33.33 %
## 33.43 %
## 33.53 %
## 33.63 %
## 33.72 %
## 33.82 %
## 33.92 %
## 34.02 %
## 34.12 %
## 34.21 %
## 34.31 %
## 34.41 %
## 34.51 %
## 34.6 %
## 34.7 %
## 34.8 %
## 34.9 %
## 35 %
## 35.09 %
## 35.19 %
## 35.29 %
## 35.39 %
## 35.48 %
## 35.58 %
## 35.68 %
## 35.78 %
## 35.87 %
## 35.97 %
## 36.07 %
## 36.17 %
## 36.27 %
## 36.36 %
## 36.46 %
## 36.56 %
## 36.66 %
## 36.75 %
## 36.85 %
## 36.95 %
## 37.05 %
## 37.15 %
## 37.24 %
## 37.34 %
## 37.44 %
## 37.54 %
## 37.63 %
## 37.73 %
## 37.83 %
## 37.93 %
## 38.03 %
## 38.12 %
## 38.22 %
## 38.32 %
## 38.42 %
## 38.51 %
## 38.61 %
## 38.71 %
## 38.81 %
## 38.91 %
## 39 %
## 39.1 %
## 39.2 %
## 39.3 %
## 39.39 %
## 39.49 %
## 39.59 %
## 39.69 %
## 39.78 %
## 39.88 %
## 39.98 %
## 40.08 %
## 40.18 %
## 40.27 %
## 40.37 %
## 40.47 %
## 40.57 %
## 40.66 %
## 40.76 %
## 40.86 %
## 40.96 %
## 41.06 %
## 41.15 %
## 41.25 %
## 41.35 %
## 41.45 %
## 41.54 %
## 41.64 %
## 41.74 %
## 41.84 %
## 41.94 %
## 42.03 %
## 42.13 %
## 42.23 %
## 42.33 %
## 42.42 %
## 42.52 %
## 42.62 %
## 42.72 %
## 42.82 %
## 42.91 %
## 43.01 %
## 43.11 %
## 43.21 %
## 43.3 %
## 43.4 %
## 43.5 %
## 43.6 %
## 43.7 %
## 43.79 %
## 43.89 %
## 43.99 %
## 44.09 %
## 44.18 %
## 44.28 %
## 44.38 %
## 44.48 %
## 44.57 %
## 44.67 %
## 44.77 %
## 44.87 %
## 44.97 %
## 45.06 %
## 45.16 %
## 45.26 %
## 45.36 %
## 45.45 %
## 45.55 %
## 45.65 %
## 45.75 %
## 45.85 %
## 45.94 %
## 46.04 %
## 46.14 %
## 46.24 %
## 46.33 %
## 46.43 %
## 46.53 %
## 46.63 %
## 46.73 %
## 46.82 %
## 46.92 %
## 47.02 %
## 47.12 %
## 47.21 %
## 47.31 %
## 47.41 %
## 47.51 %
## 47.61 %
## 47.7 %
## 47.8 %
## 47.9 %
## 48 %
## 48.09 %
## 48.19 %
## 48.29 %
## 48.39 %
## 48.48 %
## 48.58 %
## 48.68 %
## 48.78 %
## 48.88 %
## 48.97 %
## 49.07 %
## 49.17 %
## 49.27 %
## 49.36 %
## 49.46 %
## 49.56 %
## 49.66 %
## 49.76 %
## 49.85 %
## 49.95 %
## 50.05 %
## 50.15 %
## 50.24 %
## 50.34 %
## 50.44 %
## 50.54 %
## 50.64 %
## 50.73 %
## 50.83 %
## 50.93 %
## 51.03 %
## 51.12 %
## 51.22 %
## 51.32 %
## 51.42 %
## 51.52 %
## 51.61 %
## 51.71 %
## 51.81 %
## 51.91 %
## 52 %
## 52.1 %
## 52.2 %
## 52.3 %
## 52.39 %
## 52.49 %
## 52.59 %
## 52.69 %
## 52.79 %
## 52.88 %
## 52.98 %
## 53.08 %
## 53.18 %
## 53.27 %
## 53.37 %
## 53.47 %
## 53.57 %
## 53.67 %
## 53.76 %
## 53.86 %
## 53.96 %
## 54.06 %
## 54.15 %
## 54.25 %
## 54.35 %
## 54.45 %
## 54.55 %
## 54.64 %
## 54.74 %
## 54.84 %
## 54.94 %
## 55.03 %
## 55.13 %
## 55.23 %
## 55.33 %
## 55.43 %
## 55.52 %
## 55.62 %
## 55.72 %
## 55.82 %
## 55.91 %
## 56.01 %
## 56.11 %
## 56.21 %
## 56.3 %
## 56.4 %
## 56.5 %
## 56.6 %
## 56.7 %
## 56.79 %
## 56.89 %
## 56.99 %
## 57.09 %
## 57.18 %
## 57.28 %
## 57.38 %
## 57.48 %
## 57.58 %
## 57.67 %
## 57.77 %
## 57.87 %
## 57.97 %
## 58.06 %
## 58.16 %
## 58.26 %
## 58.36 %
## 58.46 %
## 58.55 %
## 58.65 %
## 58.75 %
## 58.85 %
## 58.94 %
## 59.04 %
## 59.14 %
## 59.24 %
## 59.34 %
## 59.43 %
## 59.53 %
## 59.63 %
## 59.73 %
## 59.82 %
## 59.92 %
## 60.02 %
## 60.12 %
## 60.22 %
## 60.31 %
## 60.41 %
## 60.51 %
## 60.61 %
## 60.7 %
## 60.8 %
## 60.9 %
## 61 %
## 61.09 %
## 61.19 %
## 61.29 %
## 61.39 %
## 61.49 %
## 61.58 %
## 61.68 %
## 61.78 %
## 61.88 %
## 61.97 %
## 62.07 %
## 62.17 %
## 62.27 %
## 62.37 %
## 62.46 %
## 62.56 %
## 62.66 %
## 62.76 %
## 62.85 %
## 62.95 %
## 63.05 %
## 63.15 %
## 63.25 %
## 63.34 %
## 63.44 %
## 63.54 %
## 63.64 %
## 63.73 %
## 63.83 %
## 63.93 %
## 64.03 %
## 64.13 %
## 64.22 %
## 64.32 %
## 64.42 %
## 64.52 %
## 64.61 %
## 64.71 %
## 64.81 %
## 64.91 %
## 65 %
## 65.1 %
## 65.2 %
## 65.3 %
## 65.4 %
## 65.49 %
## 65.59 %
## 65.69 %
## 65.79 %
## 65.88 %
## 65.98 %
## 66.08 %
## 66.18 %
## 66.28 %
## 66.37 %
## 66.47 %
## 66.57 %
## 66.67 %
## 66.76 %
## 66.86 %
## 66.96 %
## 67.06 %
## 67.16 %
## 67.25 %
## 67.35 %
## 67.45 %
## 67.55 %
## 67.64 %
## 67.74 %
## 67.84 %
## 67.94 %
## 68.04 %
## 68.13 %
## 68.23 %
## 68.33 %
## 68.43 %
## 68.52 %
## 68.62 %
## 68.72 %
## 68.82 %
## 68.91 %
## 69.01 %
## 69.11 %
## 69.21 %
## 69.31 %
## 69.4 %
## 69.5 %
## 69.6 %
## 69.7 %
## 69.79 %
## 69.89 %
## 69.99 %
## 70.09 %
## 70.19 %
## 70.28 %
## 70.38 %
## 70.48 %
## 70.58 %
## 70.67 %
## 70.77 %
## 70.87 %
## 70.97 %
## 71.07 %
## 71.16 %
## 71.26 %
## 71.36 %
## 71.46 %
## 71.55 %
## 71.65 %
## 71.75 %
## 71.85 %
## 71.95 %
## 72.04 %
## 72.14 %
## 72.24 %
## 72.34 %
## 72.43 %
## 72.53 %
## 72.63 %
## 72.73 %
## 72.83 %
## 72.92 %
## 73.02 %
## 73.12 %
## 73.22 %
## 73.31 %
## 73.41 %
## 73.51 %
## 73.61 %
## 73.7 %
## 73.8 %
## 73.9 %
## 74 %
## 74.1 %
## 74.19 %
## 74.29 %
## 74.39 %
## 74.49 %
## 74.58 %
## 74.68 %
## 74.78 %
## 74.88 %
## 74.98 %
## 75.07 %
## 75.17 %
## 75.27 %
## 75.37 %
## 75.46 %
## 75.56 %
## 75.66 %
## 75.76 %
## 75.86 %
## 75.95 %
## 76.05 %
## 76.15 %
## 76.25 %
## 76.34 %
## 76.44 %
## 76.54 %
## 76.64 %
## 76.74 %
## 76.83 %
## 76.93 %
## 77.03 %
## 77.13 %
## 77.22 %
## 77.32 %
## 77.42 %
## 77.52 %
## 77.61 %
## 77.71 %
## 77.81 %
## 77.91 %
## 78.01 %
## 78.1 %
## 78.2 %
## 78.3 %
## 78.4 %
## 78.49 %
## 78.59 %
## 78.69 %
## 78.79 %
## 78.89 %
## 78.98 %
## 79.08 %
## 79.18 %
## 79.28 %
## 79.37 %
## 79.47 %
## 79.57 %
## 79.67 %
## 79.77 %
## 79.86 %
## 79.96 %
## 80.06 %
## 80.16 %
## 80.25 %
## 80.35 %
## 80.45 %
## 80.55 %
## 80.65 %
## 80.74 %
## 80.84 %
## 80.94 %
## 81.04 %
## 81.13 %
## 81.23 %
## 81.33 %
## 81.43 %
## 81.52 %
## 81.62 %
## 81.72 %
## 81.82 %
## 81.92 %
## 82.01 %
## 82.11 %
## 82.21 %
## 82.31 %
## 82.4 %
## 82.5 %
## 82.6 %
## 82.7 %
## 82.8 %
## 82.89 %
## 82.99 %
## 83.09 %
## 83.19 %
## 83.28 %
## 83.38 %
## 83.48 %
## 83.58 %
## 83.68 %
## 83.77 %
## 83.87 %
## 83.97 %
## 84.07 %
## 84.16 %
## 84.26 %
## 84.36 %
## 84.46 %
## 84.56 %
## 84.65 %
## 84.75 %
## 84.85 %
## 84.95 %
## 85.04 %
## 85.14 %
## 85.24 %
## 85.34 %
## 85.43 %
## 85.53 %
## 85.63 %
## 85.73 %
## 85.83 %
## 85.92 %
## 86.02 %
## 86.12 %
## 86.22 %
## 86.31 %
## 86.41 %
## 86.51 %
## 86.61 %
## 86.71 %
## 86.8 %
## 86.9 %
## 87 %
## 87.1 %
## 87.19 %
## 87.29 %
## 87.39 %
## 87.49 %
## 87.59 %
## 87.68 %
## 87.78 %
## 87.88 %
## 87.98 %
## 88.07 %
## 88.17 %
## 88.27 %
## 88.37 %
## 88.47 %
## 88.56 %
## 88.66 %
## 88.76 %
## 88.86 %
## 88.95 %
## 89.05 %
## 89.15 %
## 89.25 %
## 89.35 %
## 89.44 %
## 89.54 %
## 89.64 %
## 89.74 %
## 89.83 %
## 89.93 %
## 90.03 %
## 90.13 %
## 90.22 %
## 90.32 %
## 90.42 %
## 90.52 %
## 90.62 %
## 90.71 %
## 90.81 %
## 90.91 %
## 91.01 %
## 91.1 %
## 91.2 %
## 91.3 %
## 91.4 %
## 91.5 %
## 91.59 %
## 91.69 %
## 91.79 %
## 91.89 %
## 91.98 %
## 92.08 %
## 92.18 %
## 92.28 %
## 92.38 %
## 92.47 %
## 92.57 %
## 92.67 %
## 92.77 %
## 92.86 %
## 92.96 %
## 93.06 %
## 93.16 %
## 93.26 %
## 93.35 %
## 93.45 %
## 93.55 %
## 93.65 %
## 93.74 %
## 93.84 %
## 93.94 %
## 94.04 %
## 94.13 %
## 94.23 %
## 94.33 %
## 94.43 %
## 94.53 %
## 94.62 %
## 94.72 %
## 94.82 %
## 94.92 %
## 95.01 %
## 95.11 %
## 95.21 %
## 95.31 %
## 95.41 %
## 95.5 %
## 95.6 %
## 95.7 %
## 95.8 %
## 95.89 %
## 95.99 %
## 96.09 %
## 96.19 %
## 96.29 %
## 96.38 %
## 96.48 %
## 96.58 %
## 96.68 %
## 96.77 %
## 96.87 %
## 96.97 %
## 97.07 %
## 97.17 %
## 97.26 %
## 97.36 %
## 97.46 %
## 97.56 %
## 97.65 %
## 97.75 %
## 97.85 %
## 97.95 %
## 98.04 %
## 98.14 %
## 98.24 %
## 98.34 %
## 98.44 %
## 98.53 %
## 98.63 %
## 98.73 %
## 98.83 %
## 98.92 %
## 99.02 %
## 99.12 %
## 99.22 %
## 99.32 %
## 99.41 %
## 99.51 %
## 99.61 %
## 99.71 %
## 99.8 %
## 99.9 %
## 100 %
## Done.
## Identifying best models (est)...
## 0.1 %
## 0.2 %
## 0.29 %
## 0.39 %
## 0.49 %
## 0.59 %
## 0.68 %
## 0.78 %
## 0.88 %
## 0.98 %
## 1.08 %
## 1.17 %
## 1.27 %
## 1.37 %
## 1.47 %
## 1.56 %
## 1.66 %
## 1.76 %
## 1.86 %
## 1.96 %
## 2.05 %
## 2.15 %
## 2.25 %
## 2.35 %
## 2.44 %
## 2.54 %
## 2.64 %
## 2.74 %
## 2.83 %
## 2.93 %
## 3.03 %
## 3.13 %
## 3.23 %
## 3.32 %
## 3.42 %
## 3.52 %
## 3.62 %
## 3.71 %
## 3.81 %
## 3.91 %
## 4.01 %
## 4.11 %
## 4.2 %
## 4.3 %
## 4.4 %
## 4.5 %
## 4.59 %
## 4.69 %
## 4.79 %
## 4.89 %
## 4.99 %
## 5.08 %
## 5.18 %
## 5.28 %
## 5.38 %
## 5.47 %
## 5.57 %
## 5.67 %
## 5.77 %
## 5.87 %
## 5.96 %
## 6.06 %
## 6.16 %
## 6.26 %
## 6.35 %
## 6.45 %
## 6.55 %
## 6.65 %
## 6.74 %
## 6.84 %
## 6.94 %
## 7.04 %
## 7.14 %
## 7.23 %
## 7.33 %
## 7.43 %
## 7.53 %
## 7.62 %
## 7.72 %
## 7.82 %
## 7.92 %
## 8.02 %
## 8.11 %
## 8.21 %
## 8.31 %
## 8.41 %
## 8.5 %
## 8.6 %
## 8.7 %
## 8.8 %
## 8.9 %
## 8.99 %
## 9.09 %
## 9.19 %
## 9.29 %
## 9.38 %
## 9.48 %
## 9.58 %
## 9.68 %
## 9.78 %
## 9.87 %
## 9.97 %
## 10.07 %
## 10.17 %
## 10.26 %
## 10.36 %
## 10.46 %
## 10.56 %
## 10.65 %
## 10.75 %
## 10.85 %
## 10.95 %
## 11.05 %
## 11.14 %
## 11.24 %
## 11.34 %
## 11.44 %
## 11.53 %
## 11.63 %
## 11.73 %
## 11.83 %
## 11.93 %
## 12.02 %
## 12.12 %
## 12.22 %
## 12.32 %
## 12.41 %
## 12.51 %
## 12.61 %
## 12.71 %
## 12.81 %
## 12.9 %
## 13 %
## 13.1 %
## 13.2 %
## 13.29 %
## 13.39 %
## 13.49 %
## 13.59 %
## 13.69 %
## 13.78 %
## 13.88 %
## 13.98 %
## 14.08 %
## 14.17 %
## 14.27 %
## 14.37 %
## 14.47 %
## 14.57 %
## 14.66 %
## 14.76 %
## 14.86 %
## 14.96 %
## 15.05 %
## 15.15 %
## 15.25 %
## 15.35 %
## 15.44 %
## 15.54 %
## 15.64 %
## 15.74 %
## 15.84 %
## 15.93 %
## 16.03 %
## 16.13 %
## 16.23 %
## 16.32 %
## 16.42 %
## 16.52 %
## 16.62 %
## 16.72 %
## 16.81 %
## 16.91 %
## 17.01 %
## 17.11 %
## 17.2 %
## 17.3 %
## 17.4 %
## 17.5 %
## 17.6 %
## 17.69 %
## 17.79 %
## 17.89 %
## 17.99 %
## 18.08 %
## 18.18 %
## 18.28 %
## 18.38 %
## 18.48 %
## 18.57 %
## 18.67 %
## 18.77 %
## 18.87 %
## 18.96 %
## 19.06 %
## 19.16 %
## 19.26 %
## 19.35 %
## 19.45 %
## 19.55 %
## 19.65 %
## 19.75 %
## 19.84 %
## 19.94 %
## 20.04 %
## 20.14 %
## 20.23 %
## 20.33 %
## 20.43 %
## 20.53 %
## 20.63 %
## 20.72 %
## 20.82 %
## 20.92 %
## 21.02 %
## 21.11 %
## 21.21 %
## 21.31 %
## 21.41 %
## 21.51 %
## 21.6 %
## 21.7 %
## 21.8 %
## 21.9 %
## 21.99 %
## 22.09 %
## 22.19 %
## 22.29 %
## 22.39 %
## 22.48 %
## 22.58 %
## 22.68 %
## 22.78 %
## 22.87 %
## 22.97 %
## 23.07 %
## 23.17 %
## 23.26 %
## 23.36 %
## 23.46 %
## 23.56 %
## 23.66 %
## 23.75 %
## 23.85 %
## 23.95 %
## 24.05 %
## 24.14 %
## 24.24 %
## 24.34 %
## 24.44 %
## 24.54 %
## 24.63 %
## 24.73 %
## 24.83 %
## 24.93 %
## 25.02 %
## 25.12 %
## 25.22 %
## 25.32 %
## 25.42 %
## 25.51 %
## 25.61 %
## 25.71 %
## 25.81 %
## 25.9 %
## 26 %
## 26.1 %
## 26.2 %
## 26.3 %
## 26.39 %
## 26.49 %
## 26.59 %
## 26.69 %
## 26.78 %
## 26.88 %
## 26.98 %
## 27.08 %
## 27.17 %
## 27.27 %
## 27.37 %
## 27.47 %
## 27.57 %
## 27.66 %
## 27.76 %
## 27.86 %
## 27.96 %
## 28.05 %
## 28.15 %
## 28.25 %
## 28.35 %
## 28.45 %
## 28.54 %
## 28.64 %
## 28.74 %
## 28.84 %
## 28.93 %
## 29.03 %
## 29.13 %
## 29.23 %
## 29.33 %
## 29.42 %
## 29.52 %
## 29.62 %
## 29.72 %
## 29.81 %
## 29.91 %
## 30.01 %
## 30.11 %
## 30.21 %
## 30.3 %
## 30.4 %
## 30.5 %
## 30.6 %
## 30.69 %
## 30.79 %
## 30.89 %
## 30.99 %
## 31.09 %
## 31.18 %
## 31.28 %
## 31.38 %
## 31.48 %
## 31.57 %
## 31.67 %
## 31.77 %
## 31.87 %
## 31.96 %
## 32.06 %
## 32.16 %
## 32.26 %
## 32.36 %
## 32.45 %
## 32.55 %
## 32.65 %
## 32.75 %
## 32.84 %
## 32.94 %
## 33.04 %
## 33.14 %
## 33.24 %
## 33.33 %
## 33.43 %
## 33.53 %
## 33.63 %
## 33.72 %
## 33.82 %
## 33.92 %
## 34.02 %
## 34.12 %
## 34.21 %
## 34.31 %
## 34.41 %
## 34.51 %
## 34.6 %
## 34.7 %
## 34.8 %
## 34.9 %
## 35 %
## 35.09 %
## 35.19 %
## 35.29 %
## 35.39 %
## 35.48 %
## 35.58 %
## 35.68 %
## 35.78 %
## 35.87 %
## 35.97 %
## 36.07 %
## 36.17 %
## 36.27 %
## 36.36 %
## 36.46 %
## 36.56 %
## 36.66 %
## 36.75 %
## 36.85 %
## 36.95 %
## 37.05 %
## 37.15 %
## 37.24 %
## 37.34 %
## 37.44 %
## 37.54 %
## 37.63 %
## 37.73 %
## 37.83 %
## 37.93 %
## 38.03 %
## 38.12 %
## 38.22 %
## 38.32 %
## 38.42 %
## 38.51 %
## 38.61 %
## 38.71 %
## 38.81 %
## 38.91 %
## 39 %
## 39.1 %
## 39.2 %
## 39.3 %
## 39.39 %
## 39.49 %
## 39.59 %
## 39.69 %
## 39.78 %
## 39.88 %
## 39.98 %
## 40.08 %
## 40.18 %
## 40.27 %
## 40.37 %
## 40.47 %
## 40.57 %
## 40.66 %
## 40.76 %
## 40.86 %
## 40.96 %
## 41.06 %
## 41.15 %
## 41.25 %
## 41.35 %
## 41.45 %
## 41.54 %
## 41.64 %
## 41.74 %
## 41.84 %
## 41.94 %
## 42.03 %
## 42.13 %
## 42.23 %
## 42.33 %
## 42.42 %
## 42.52 %
## 42.62 %
## 42.72 %
## 42.82 %
## 42.91 %
## 43.01 %
## 43.11 %
## 43.21 %
## 43.3 %
## 43.4 %
## 43.5 %
## 43.6 %
## 43.7 %
## 43.79 %
## 43.89 %
## 43.99 %
## 44.09 %
## 44.18 %
## 44.28 %
## 44.38 %
## 44.48 %
## 44.57 %
## 44.67 %
## 44.77 %
## 44.87 %
## 44.97 %
## 45.06 %
## 45.16 %
## 45.26 %
## 45.36 %
## 45.45 %
## 45.55 %
## 45.65 %
## 45.75 %
## 45.85 %
## 45.94 %
## 46.04 %
## 46.14 %
## 46.24 %
## 46.33 %
## 46.43 %
## 46.53 %
## 46.63 %
## 46.73 %
## 46.82 %
## 46.92 %
## 47.02 %
## 47.12 %
## 47.21 %
## 47.31 %
## 47.41 %
## 47.51 %
## 47.61 %
## 47.7 %
## 47.8 %
## 47.9 %
## 48 %
## 48.09 %
## 48.19 %
## 48.29 %
## 48.39 %
## 48.48 %
## 48.58 %
## 48.68 %
## 48.78 %
## 48.88 %
## 48.97 %
## 49.07 %
## 49.17 %
## 49.27 %
## 49.36 %
## 49.46 %
## 49.56 %
## 49.66 %
## 49.76 %
## 49.85 %
## 49.95 %
## 50.05 %
## 50.15 %
## 50.24 %
## 50.34 %
## 50.44 %
## 50.54 %
## 50.64 %
## 50.73 %
## 50.83 %
## 50.93 %
## 51.03 %
## 51.12 %
## 51.22 %
## 51.32 %
## 51.42 %
## 51.52 %
## 51.61 %
## 51.71 %
## 51.81 %
## 51.91 %
## 52 %
## 52.1 %
## 52.2 %
## 52.3 %
## 52.39 %
## 52.49 %
## 52.59 %
## 52.69 %
## 52.79 %
## 52.88 %
## 52.98 %
## 53.08 %
## 53.18 %
## 53.27 %
## 53.37 %
## 53.47 %
## 53.57 %
## 53.67 %
## 53.76 %
## 53.86 %
## 53.96 %
## 54.06 %
## 54.15 %
## 54.25 %
## 54.35 %
## 54.45 %
## 54.55 %
## 54.64 %
## 54.74 %
## 54.84 %
## 54.94 %
## 55.03 %
## 55.13 %
## 55.23 %
## 55.33 %
## 55.43 %
## 55.52 %
## 55.62 %
## 55.72 %
## 55.82 %
## 55.91 %
## 56.01 %
## 56.11 %
## 56.21 %
## 56.3 %
## 56.4 %
## 56.5 %
## 56.6 %
## 56.7 %
## 56.79 %
## 56.89 %
## 56.99 %
## 57.09 %
## 57.18 %
## 57.28 %
## 57.38 %
## 57.48 %
## 57.58 %
## 57.67 %
## 57.77 %
## 57.87 %
## 57.97 %
## 58.06 %
## 58.16 %
## 58.26 %
## 58.36 %
## 58.46 %
## 58.55 %
## 58.65 %
## 58.75 %
## 58.85 %
## 58.94 %
## 59.04 %
## 59.14 %
## 59.24 %
## 59.34 %
## 59.43 %
## 59.53 %
## 59.63 %
## 59.73 %
## 59.82 %
## 59.92 %
## 60.02 %
## 60.12 %
## 60.22 %
## 60.31 %
## 60.41 %
## 60.51 %
## 60.61 %
## 60.7 %
## 60.8 %
## 60.9 %
## 61 %
## 61.09 %
## 61.19 %
## 61.29 %
## 61.39 %
## 61.49 %
## 61.58 %
## 61.68 %
## 61.78 %
## 61.88 %
## 61.97 %
## 62.07 %
## 62.17 %
## 62.27 %
## 62.37 %
## 62.46 %
## 62.56 %
## 62.66 %
## 62.76 %
## 62.85 %
## 62.95 %
## 63.05 %
## 63.15 %
## 63.25 %
## 63.34 %
## 63.44 %
## 63.54 %
## 63.64 %
## 63.73 %
## 63.83 %
## 63.93 %
## 64.03 %
## 64.13 %
## 64.22 %
## 64.32 %
## 64.42 %
## 64.52 %
## 64.61 %
## 64.71 %
## 64.81 %
## 64.91 %
## 65 %
## 65.1 %
## 65.2 %
## 65.3 %
## 65.4 %
## 65.49 %
## 65.59 %
## 65.69 %
## 65.79 %
## 65.88 %
## 65.98 %
## 66.08 %
## 66.18 %
## 66.28 %
## 66.37 %
## 66.47 %
## 66.57 %
## 66.67 %
## 66.76 %
## 66.86 %
## 66.96 %
## 67.06 %
## 67.16 %
## 67.25 %
## 67.35 %
## 67.45 %
## 67.55 %
## 67.64 %
## 67.74 %
## 67.84 %
## 67.94 %
## 68.04 %
## 68.13 %
## 68.23 %
## 68.33 %
## 68.43 %
## 68.52 %
## 68.62 %
## 68.72 %
## 68.82 %
## 68.91 %
## 69.01 %
## 69.11 %
## 69.21 %
## 69.31 %
## 69.4 %
## 69.5 %
## 69.6 %
## 69.7 %
## 69.79 %
## 69.89 %
## 69.99 %
## 70.09 %
## 70.19 %
## 70.28 %
## 70.38 %
## 70.48 %
## 70.58 %
## 70.67 %
## 70.77 %
## 70.87 %
## 70.97 %
## 71.07 %
## 71.16 %
## 71.26 %
## 71.36 %
## 71.46 %
## 71.55 %
## 71.65 %
## 71.75 %
## 71.85 %
## 71.95 %
## 72.04 %
## 72.14 %
## 72.24 %
## 72.34 %
## 72.43 %
## 72.53 %
## 72.63 %
## 72.73 %
## 72.83 %
## 72.92 %
## 73.02 %
## 73.12 %
## 73.22 %
## 73.31 %
## 73.41 %
## 73.51 %
## 73.61 %
## 73.7 %
## 73.8 %
## 73.9 %
## 74 %
## 74.1 %
## 74.19 %
## 74.29 %
## 74.39 %
## 74.49 %
## 74.58 %
## 74.68 %
## 74.78 %
## 74.88 %
## 74.98 %
## 75.07 %
## 75.17 %
## 75.27 %
## 75.37 %
## 75.46 %
## 75.56 %
## 75.66 %
## 75.76 %
## 75.86 %
## 75.95 %
## 76.05 %
## 76.15 %
## 76.25 %
## 76.34 %
## 76.44 %
## 76.54 %
## 76.64 %
## 76.74 %
## 76.83 %
## 76.93 %
## 77.03 %
## 77.13 %
## 77.22 %
## 77.32 %
## 77.42 %
## 77.52 %
## 77.61 %
## 77.71 %
## 77.81 %
## 77.91 %
## 78.01 %
## 78.1 %
## 78.2 %
## 78.3 %
## 78.4 %
## 78.49 %
## 78.59 %
## 78.69 %
## 78.79 %
## 78.89 %
## 78.98 %
## 79.08 %
## 79.18 %
## 79.28 %
## 79.37 %
## 79.47 %
## 79.57 %
## 79.67 %
## 79.77 %
## 79.86 %
## 79.96 %
## 80.06 %
## 80.16 %
## 80.25 %
## 80.35 %
## 80.45 %
## 80.55 %
## 80.65 %
## 80.74 %
## 80.84 %
## 80.94 %
## 81.04 %
## 81.13 %
## 81.23 %
## 81.33 %
## 81.43 %
## 81.52 %
## 81.62 %
## 81.72 %
## 81.82 %
## 81.92 %
## 82.01 %
## 82.11 %
## 82.21 %
## 82.31 %
## 82.4 %
## 82.5 %
## 82.6 %
## 82.7 %
## 82.8 %
## 82.89 %
## 82.99 %
## 83.09 %
## 83.19 %
## 83.28 %
## 83.38 %
## 83.48 %
## 83.58 %
## 83.68 %
## 83.77 %
## 83.87 %
## 83.97 %
## 84.07 %
## 84.16 %
## 84.26 %
## 84.36 %
## 84.46 %
## 84.56 %
## 84.65 %
## 84.75 %
## 84.85 %
## 84.95 %
## 85.04 %
## 85.14 %
## 85.24 %
## 85.34 %
## 85.43 %
## 85.53 %
## 85.63 %
## 85.73 %
## 85.83 %
## 85.92 %
## 86.02 %
## 86.12 %
## 86.22 %
## 86.31 %
## 86.41 %
## 86.51 %
## 86.61 %
## 86.71 %
## 86.8 %
## 86.9 %
## 87 %
## 87.1 %
## 87.19 %
## 87.29 %
## 87.39 %
## 87.49 %
## 87.59 %
## 87.68 %
## 87.78 %
## 87.88 %
## 87.98 %
## 88.07 %
## 88.17 %
## 88.27 %
## 88.37 %
## 88.47 %
## 88.56 %
## 88.66 %
## 88.76 %
## 88.86 %
## 88.95 %
## 89.05 %
## 89.15 %
## 89.25 %
## 89.35 %
## 89.44 %
## 89.54 %
## 89.64 %
## 89.74 %
## 89.83 %
## 89.93 %
## 90.03 %
## 90.13 %
## 90.22 %
## 90.32 %
## 90.42 %
## 90.52 %
## 90.62 %
## 90.71 %
## 90.81 %
## 90.91 %
## 91.01 %
## 91.1 %
## 91.2 %
## 91.3 %
## 91.4 %
## 91.5 %
## 91.59 %
## 91.69 %
## 91.79 %
## 91.89 %
## 91.98 %
## 92.08 %
## 92.18 %
## 92.28 %
## 92.38 %
## 92.47 %
## 92.57 %
## 92.67 %
## 92.77 %
## 92.86 %
## 92.96 %
## 93.06 %
## 93.16 %
## 93.26 %
## 93.35 %
## 93.45 %
## 93.55 %
## 93.65 %
## 93.74 %
## 93.84 %
## 93.94 %
## 94.04 %
## 94.13 %
## 94.23 %
## 94.33 %
## 94.43 %
## 94.53 %
## 94.62 %
## 94.72 %
## 94.82 %
## 94.92 %
## 95.01 %
## 95.11 %
## 95.21 %
## 95.31 %
## 95.41 %
## 95.5 %
## 95.6 %
## 95.7 %
## 95.8 %
## 95.89 %
## 95.99 %
## 96.09 %
## 96.19 %
## 96.29 %
## 96.38 %
## 96.48 %
## 96.58 %
## 96.68 %
## 96.77 %
## 96.87 %
## 96.97 %
## 97.07 %
## 97.17 %
## 97.26 %
## 97.36 %
## 97.46 %
## 97.56 %
## 97.65 %
## 97.75 %
## 97.85 %
## 97.95 %
## 98.04 %
## 98.14 %
## 98.24 %
## 98.34 %
## 98.44 %
## 98.53 %
## 98.63 %
## 98.73 %
## 98.83 %
## 98.92 %
## 99.02 %
## 99.12 %
## 99.22 %
## 99.32 %
## 99.41 %
## 99.51 %
## 99.61 %
## 99.71 %
## 99.8 %
## 99.9 %
## 100 %
## Binding data...
## Arranging data...
## Identifying best models (p)...
## 0.1 %
## 0.2 %
## 0.29 %
## 0.39 %
## 0.49 %
## 0.59 %
## 0.68 %
## 0.78 %
## 0.88 %
## 0.98 %
## 1.08 %
## 1.17 %
## 1.27 %
## 1.37 %
## 1.47 %
## 1.56 %
## 1.66 %
## 1.76 %
## 1.86 %
## 1.96 %
## 2.05 %
## 2.15 %
## 2.25 %
## 2.35 %
## 2.44 %
## 2.54 %
## 2.64 %
## 2.74 %
## 2.83 %
## 2.93 %
## 3.03 %
## 3.13 %
## 3.23 %
## 3.32 %
## 3.42 %
## 3.52 %
## 3.62 %
## 3.71 %
## 3.81 %
## 3.91 %
## 4.01 %
## 4.11 %
## 4.2 %
## 4.3 %
## 4.4 %
## 4.5 %
## 4.59 %
## 4.69 %
## 4.79 %
## 4.89 %
## 4.99 %
## 5.08 %
## 5.18 %
## 5.28 %
## 5.38 %
## 5.47 %
## 5.57 %
## 5.67 %
## 5.77 %
## 5.87 %
## 5.96 %
## 6.06 %
## 6.16 %
## 6.26 %
## 6.35 %
## 6.45 %
## 6.55 %
## 6.65 %
## 6.74 %
## 6.84 %
## 6.94 %
## 7.04 %
## 7.14 %
## 7.23 %
## 7.33 %
## 7.43 %
## 7.53 %
## 7.62 %
## 7.72 %
## 7.82 %
## 7.92 %
## 8.02 %
## 8.11 %
## 8.21 %
## 8.31 %
## 8.41 %
## 8.5 %
## 8.6 %
## 8.7 %
## 8.8 %
## 8.9 %
## 8.99 %
## 9.09 %
## 9.19 %
## 9.29 %
## 9.38 %
## 9.48 %
## 9.58 %
## 9.68 %
## 9.78 %
## 9.87 %
## 9.97 %
## 10.07 %
## 10.17 %
## 10.26 %
## 10.36 %
## 10.46 %
## 10.56 %
## 10.65 %
## 10.75 %
## 10.85 %
## 10.95 %
## 11.05 %
## 11.14 %
## 11.24 %
## 11.34 %
## 11.44 %
## 11.53 %
## 11.63 %
## 11.73 %
## 11.83 %
## 11.93 %
## 12.02 %
## 12.12 %
## 12.22 %
## 12.32 %
## 12.41 %
## 12.51 %
## 12.61 %
## 12.71 %
## 12.81 %
## 12.9 %
## 13 %
## 13.1 %
## 13.2 %
## 13.29 %
## 13.39 %
## 13.49 %
## 13.59 %
## 13.69 %
## 13.78 %
## 13.88 %
## 13.98 %
## 14.08 %
## 14.17 %
## 14.27 %
## 14.37 %
## 14.47 %
## 14.57 %
## 14.66 %
## 14.76 %
## 14.86 %
## 14.96 %
## 15.05 %
## 15.15 %
## 15.25 %
## 15.35 %
## 15.44 %
## 15.54 %
## 15.64 %
## 15.74 %
## 15.84 %
## 15.93 %
## 16.03 %
## 16.13 %
## 16.23 %
## 16.32 %
## 16.42 %
## 16.52 %
## 16.62 %
## 16.72 %
## 16.81 %
## 16.91 %
## 17.01 %
## 17.11 %
## 17.2 %
## 17.3 %
## 17.4 %
## 17.5 %
## 17.6 %
## 17.69 %
## 17.79 %
## 17.89 %
## 17.99 %
## 18.08 %
## 18.18 %
## 18.28 %
## 18.38 %
## 18.48 %
## 18.57 %
## 18.67 %
## 18.77 %
## 18.87 %
## 18.96 %
## 19.06 %
## 19.16 %
## 19.26 %
## 19.35 %
## 19.45 %
## 19.55 %
## 19.65 %
## 19.75 %
## 19.84 %
## 19.94 %
## 20.04 %
## 20.14 %
## 20.23 %
## 20.33 %
## 20.43 %
## 20.53 %
## 20.63 %
## 20.72 %
## 20.82 %
## 20.92 %
## 21.02 %
## 21.11 %
## 21.21 %
## 21.31 %
## 21.41 %
## 21.51 %
## 21.6 %
## 21.7 %
## 21.8 %
## 21.9 %
## 21.99 %
## 22.09 %
## 22.19 %
## 22.29 %
## 22.39 %
## 22.48 %
## 22.58 %
## 22.68 %
## 22.78 %
## 22.87 %
## 22.97 %
## 23.07 %
## 23.17 %
## 23.26 %
## 23.36 %
## 23.46 %
## 23.56 %
## 23.66 %
## 23.75 %
## 23.85 %
## 23.95 %
## 24.05 %
## 24.14 %
## 24.24 %
## 24.34 %
## 24.44 %
## 24.54 %
## 24.63 %
## 24.73 %
## 24.83 %
## 24.93 %
## 25.02 %
## 25.12 %
## 25.22 %
## 25.32 %
## 25.42 %
## 25.51 %
## 25.61 %
## 25.71 %
## 25.81 %
## 25.9 %
## 26 %
## 26.1 %
## 26.2 %
## 26.3 %
## 26.39 %
## 26.49 %
## 26.59 %
## 26.69 %
## 26.78 %
## 26.88 %
## 26.98 %
## 27.08 %
## 27.17 %
## 27.27 %
## 27.37 %
## 27.47 %
## 27.57 %
## 27.66 %
## 27.76 %
## 27.86 %
## 27.96 %
## 28.05 %
## 28.15 %
## 28.25 %
## 28.35 %
## 28.45 %
## 28.54 %
## 28.64 %
## 28.74 %
## 28.84 %
## 28.93 %
## 29.03 %
## 29.13 %
## 29.23 %
## 29.33 %
## 29.42 %
## 29.52 %
## 29.62 %
## 29.72 %
## 29.81 %
## 29.91 %
## 30.01 %
## 30.11 %
## 30.21 %
## 30.3 %
## 30.4 %
## 30.5 %
## 30.6 %
## 30.69 %
## 30.79 %
## 30.89 %
## 30.99 %
## 31.09 %
## 31.18 %
## 31.28 %
## 31.38 %
## 31.48 %
## 31.57 %
## 31.67 %
## 31.77 %
## 31.87 %
## 31.96 %
## 32.06 %
## 32.16 %
## 32.26 %
## 32.36 %
## 32.45 %
## 32.55 %
## 32.65 %
## 32.75 %
## 32.84 %
## 32.94 %
## 33.04 %
## 33.14 %
## 33.24 %
## 33.33 %
## 33.43 %
## 33.53 %
## 33.63 %
## 33.72 %
## 33.82 %
## 33.92 %
## 34.02 %
## 34.12 %
## 34.21 %
## 34.31 %
## 34.41 %
## 34.51 %
## 34.6 %
## 34.7 %
## 34.8 %
## 34.9 %
## 35 %
## 35.09 %
## 35.19 %
## 35.29 %
## 35.39 %
## 35.48 %
## 35.58 %
## 35.68 %
## 35.78 %
## 35.87 %
## 35.97 %
## 36.07 %
## 36.17 %
## 36.27 %
## 36.36 %
## 36.46 %
## 36.56 %
## 36.66 %
## 36.75 %
## 36.85 %
## 36.95 %
## 37.05 %
## 37.15 %
## 37.24 %
## 37.34 %
## 37.44 %
## 37.54 %
## 37.63 %
## 37.73 %
## 37.83 %
## 37.93 %
## 38.03 %
## 38.12 %
## 38.22 %
## 38.32 %
## 38.42 %
## 38.51 %
## 38.61 %
## 38.71 %
## 38.81 %
## 38.91 %
## 39 %
## 39.1 %
## 39.2 %
## 39.3 %
## 39.39 %
## 39.49 %
## 39.59 %
## 39.69 %
## 39.78 %
## 39.88 %
## 39.98 %
## 40.08 %
## 40.18 %
## 40.27 %
## 40.37 %
## 40.47 %
## 40.57 %
## 40.66 %
## 40.76 %
## 40.86 %
## 40.96 %
## 41.06 %
## 41.15 %
## 41.25 %
## 41.35 %
## 41.45 %
## 41.54 %
## 41.64 %
## 41.74 %
## 41.84 %
## 41.94 %
## 42.03 %
## 42.13 %
## 42.23 %
## 42.33 %
## 42.42 %
## 42.52 %
## 42.62 %
## 42.72 %
## 42.82 %
## 42.91 %
## 43.01 %
## 43.11 %
## 43.21 %
## 43.3 %
## 43.4 %
## 43.5 %
## 43.6 %
## 43.7 %
## 43.79 %
## 43.89 %
## 43.99 %
## 44.09 %
## 44.18 %
## 44.28 %
## 44.38 %
## 44.48 %
## 44.57 %
## 44.67 %
## 44.77 %
## 44.87 %
## 44.97 %
## 45.06 %
## 45.16 %
## 45.26 %
## 45.36 %
## 45.45 %
## 45.55 %
## 45.65 %
## 45.75 %
## 45.85 %
## 45.94 %
## 46.04 %
## 46.14 %
## 46.24 %
## 46.33 %
## 46.43 %
## 46.53 %
## 46.63 %
## 46.73 %
## 46.82 %
## 46.92 %
## 47.02 %
## 47.12 %
## 47.21 %
## 47.31 %
## 47.41 %
## 47.51 %
## 47.61 %
## 47.7 %
## 47.8 %
## 47.9 %
## 48 %
## 48.09 %
## 48.19 %
## 48.29 %
## 48.39 %
## 48.48 %
## 48.58 %
## 48.68 %
## 48.78 %
## 48.88 %
## 48.97 %
## 49.07 %
## 49.17 %
## 49.27 %
## 49.36 %
## 49.46 %
## 49.56 %
## 49.66 %
## 49.76 %
## 49.85 %
## 49.95 %
## 50.05 %
## 50.15 %
## 50.24 %
## 50.34 %
## 50.44 %
## 50.54 %
## 50.64 %
## 50.73 %
## 50.83 %
## 50.93 %
## 51.03 %
## 51.12 %
## 51.22 %
## 51.32 %
## 51.42 %
## 51.52 %
## 51.61 %
## 51.71 %
## 51.81 %
## 51.91 %
## 52 %
## 52.1 %
## 52.2 %
## 52.3 %
## 52.39 %
## 52.49 %
## 52.59 %
## 52.69 %
## 52.79 %
## 52.88 %
## 52.98 %
## 53.08 %
## 53.18 %
## 53.27 %
## 53.37 %
## 53.47 %
## 53.57 %
## 53.67 %
## 53.76 %
## 53.86 %
## 53.96 %
## 54.06 %
## 54.15 %
## 54.25 %
## 54.35 %
## 54.45 %
## 54.55 %
## 54.64 %
## 54.74 %
## 54.84 %
## 54.94 %
## 55.03 %
## 55.13 %
## 55.23 %
## 55.33 %
## 55.43 %
## 55.52 %
## 55.62 %
## 55.72 %
## 55.82 %
## 55.91 %
## 56.01 %
## 56.11 %
## 56.21 %
## 56.3 %
## 56.4 %
## 56.5 %
## 56.6 %
## 56.7 %
## 56.79 %
## 56.89 %
## 56.99 %
## 57.09 %
## 57.18 %
## 57.28 %
## 57.38 %
## 57.48 %
## 57.58 %
## 57.67 %
## 57.77 %
## 57.87 %
## 57.97 %
## 58.06 %
## 58.16 %
## 58.26 %
## 58.36 %
## 58.46 %
## 58.55 %
## 58.65 %
## 58.75 %
## 58.85 %
## 58.94 %
## 59.04 %
## 59.14 %
## 59.24 %
## 59.34 %
## 59.43 %
## 59.53 %
## 59.63 %
## 59.73 %
## 59.82 %
## 59.92 %
## 60.02 %
## 60.12 %
## 60.22 %
## 60.31 %
## 60.41 %
## 60.51 %
## 60.61 %
## 60.7 %
## 60.8 %
## 60.9 %
## 61 %
## 61.09 %
## 61.19 %
## 61.29 %
## 61.39 %
## 61.49 %
## 61.58 %
## 61.68 %
## 61.78 %
## 61.88 %
## 61.97 %
## 62.07 %
## 62.17 %
## 62.27 %
## 62.37 %
## 62.46 %
## 62.56 %
## 62.66 %
## 62.76 %
## 62.85 %
## 62.95 %
## 63.05 %
## 63.15 %
## 63.25 %
## 63.34 %
## 63.44 %
## 63.54 %
## 63.64 %
## 63.73 %
## 63.83 %
## 63.93 %
## 64.03 %
## 64.13 %
## 64.22 %
## 64.32 %
## 64.42 %
## 64.52 %
## 64.61 %
## 64.71 %
## 64.81 %
## 64.91 %
## 65 %
## 65.1 %
## 65.2 %
## 65.3 %
## 65.4 %
## 65.49 %
## 65.59 %
## 65.69 %
## 65.79 %
## 65.88 %
## 65.98 %
## 66.08 %
## 66.18 %
## 66.28 %
## 66.37 %
## 66.47 %
## 66.57 %
## 66.67 %
## 66.76 %
## 66.86 %
## 66.96 %
## 67.06 %
## 67.16 %
## 67.25 %
## 67.35 %
## 67.45 %
## 67.55 %
## 67.64 %
## 67.74 %
## 67.84 %
## 67.94 %
## 68.04 %
## 68.13 %
## 68.23 %
## 68.33 %
## 68.43 %
## 68.52 %
## 68.62 %
## 68.72 %
## 68.82 %
## 68.91 %
## 69.01 %
## 69.11 %
## 69.21 %
## 69.31 %
## 69.4 %
## 69.5 %
## 69.6 %
## 69.7 %
## 69.79 %
## 69.89 %
## 69.99 %
## 70.09 %
## 70.19 %
## 70.28 %
## 70.38 %
## 70.48 %
## 70.58 %
## 70.67 %
## 70.77 %
## 70.87 %
## 70.97 %
## 71.07 %
## 71.16 %
## 71.26 %
## 71.36 %
## 71.46 %
## 71.55 %
## 71.65 %
## 71.75 %
## 71.85 %
## 71.95 %
## 72.04 %
## 72.14 %
## 72.24 %
## 72.34 %
## 72.43 %
## 72.53 %
## 72.63 %
## 72.73 %
## 72.83 %
## 72.92 %
## 73.02 %
## 73.12 %
## 73.22 %
## 73.31 %
## 73.41 %
## 73.51 %
## 73.61 %
## 73.7 %
## 73.8 %
## 73.9 %
## 74 %
## 74.1 %
## 74.19 %
## 74.29 %
## 74.39 %
## 74.49 %
## 74.58 %
## 74.68 %
## 74.78 %
## 74.88 %
## 74.98 %
## 75.07 %
## 75.17 %
## 75.27 %
## 75.37 %
## 75.46 %
## 75.56 %
## 75.66 %
## 75.76 %
## 75.86 %
## 75.95 %
## 76.05 %
## 76.15 %
## 76.25 %
## 76.34 %
## 76.44 %
## 76.54 %
## 76.64 %
## 76.74 %
## 76.83 %
## 76.93 %
## 77.03 %
## 77.13 %
## 77.22 %
## 77.32 %
## 77.42 %
## 77.52 %
## 77.61 %
## 77.71 %
## 77.81 %
## 77.91 %
## 78.01 %
## 78.1 %
## 78.2 %
## 78.3 %
## 78.4 %
## 78.49 %
## 78.59 %
## 78.69 %
## 78.79 %
## 78.89 %
## 78.98 %
## 79.08 %
## 79.18 %
## 79.28 %
## 79.37 %
## 79.47 %
## 79.57 %
## 79.67 %
## 79.77 %
## 79.86 %
## 79.96 %
## 80.06 %
## 80.16 %
## 80.25 %
## 80.35 %
## 80.45 %
## 80.55 %
## 80.65 %
## 80.74 %
## 80.84 %
## 80.94 %
## 81.04 %
## 81.13 %
## 81.23 %
## 81.33 %
## 81.43 %
## 81.52 %
## 81.62 %
## 81.72 %
## 81.82 %
## 81.92 %
## 82.01 %
## 82.11 %
## 82.21 %
## 82.31 %
## 82.4 %
## 82.5 %
## 82.6 %
## 82.7 %
## 82.8 %
## 82.89 %
## 82.99 %
## 83.09 %
## 83.19 %
## 83.28 %
## 83.38 %
## 83.48 %
## 83.58 %
## 83.68 %
## 83.77 %
## 83.87 %
## 83.97 %
## 84.07 %
## 84.16 %
## 84.26 %
## 84.36 %
## 84.46 %
## 84.56 %
## 84.65 %
## 84.75 %
## 84.85 %
## 84.95 %
## 85.04 %
## 85.14 %
## 85.24 %
## 85.34 %
## 85.43 %
## 85.53 %
## 85.63 %
## 85.73 %
## 85.83 %
## 85.92 %
## 86.02 %
## 86.12 %
## 86.22 %
## 86.31 %
## 86.41 %
## 86.51 %
## 86.61 %
## 86.71 %
## 86.8 %
## 86.9 %
## 87 %
## 87.1 %
## 87.19 %
## 87.29 %
## 87.39 %
## 87.49 %
## 87.59 %
## 87.68 %
## 87.78 %
## 87.88 %
## 87.98 %
## 88.07 %
## 88.17 %
## 88.27 %
## 88.37 %
## 88.47 %
## 88.56 %
## 88.66 %
## 88.76 %
## 88.86 %
## 88.95 %
## 89.05 %
## 89.15 %
## 89.25 %
## 89.35 %
## 89.44 %
## 89.54 %
## 89.64 %
## 89.74 %
## 89.83 %
## 89.93 %
## 90.03 %
## 90.13 %
## 90.22 %
## 90.32 %
## 90.42 %
## 90.52 %
## 90.62 %
## 90.71 %
## 90.81 %
## 90.91 %
## 91.01 %
## 91.1 %
## 91.2 %
## 91.3 %
## 91.4 %
## 91.5 %
## 91.59 %
## 91.69 %
## 91.79 %
## 91.89 %
## 91.98 %
## 92.08 %
## 92.18 %
## 92.28 %
## 92.38 %
## 92.47 %
## 92.57 %
## 92.67 %
## 92.77 %
## 92.86 %
## 92.96 %
## 93.06 %
## 93.16 %
## 93.26 %
## 93.35 %
## 93.45 %
## 93.55 %
## 93.65 %
## 93.74 %
## 93.84 %
## 93.94 %
## 94.04 %
## 94.13 %
## 94.23 %
## 94.33 %
## 94.43 %
## 94.53 %
## 94.62 %
## 94.72 %
## 94.82 %
## 94.92 %
## 95.01 %
## 95.11 %
## 95.21 %
## 95.31 %
## 95.41 %
## 95.5 %
## 95.6 %
## 95.7 %
## 95.8 %
## 95.89 %
## 95.99 %
## 96.09 %
## 96.19 %
## 96.29 %
## 96.38 %
## 96.48 %
## 96.58 %
## 96.68 %
## 96.77 %
## 96.87 %
## 96.97 %
## 97.07 %
## 97.17 %
## 97.26 %
## 97.36 %
## 97.46 %
## 97.56 %
## 97.65 %
## 97.75 %
## 97.85 %
## 97.95 %
## 98.04 %
## 98.14 %
## 98.24 %
## 98.34 %
## 98.44 %
## 98.53 %
## 98.63 %
## 98.73 %
## 98.83 %
## 98.92 %
## 99.02 %
## 99.12 %
## 99.22 %
## 99.32 %
## 99.41 %
## 99.51 %
## 99.61 %
## 99.71 %
## 99.8 %
## 99.9 %
## 100 %
## Binding data...
## Arranging data...
## Done.
## Extracting results...
## 0.1 %
## 0.2 %
## 0.29 %
## 0.39 %
## 0.49 %
## 0.59 %
## 0.68 %
## 0.78 %
## 0.88 %
## 0.98 %
## 1.08 %
## 1.17 %
## 1.27 %
## 1.37 %
## 1.47 %
## 1.56 %
## 1.66 %
## 1.76 %
## 1.86 %
## 1.96 %
## 2.05 %
## 2.15 %
## 2.25 %
## 2.35 %
## 2.44 %
## 2.54 %
## 2.64 %
## 2.74 %
## 2.83 %
## 2.93 %
## 3.03 %
## 3.13 %
## 3.23 %
## 3.32 %
## 3.42 %
## 3.52 %
## 3.62 %
## 3.71 %
## 3.81 %
## 3.91 %
## 4.01 %
## 4.11 %
## 4.2 %
## 4.3 %
## 4.4 %
## 4.5 %
## 4.59 %
## 4.69 %
## 4.79 %
## 4.89 %
## 4.99 %
## 5.08 %
## 5.18 %
## 5.28 %
## 5.38 %
## 5.47 %
## 5.57 %
## 5.67 %
## 5.77 %
## 5.87 %
## 5.96 %
## 6.06 %
## 6.16 %
## 6.26 %
## 6.35 %
## 6.45 %
## 6.55 %
## 6.65 %
## 6.74 %
## 6.84 %
## 6.94 %
## 7.04 %
## 7.14 %
## 7.23 %
## 7.33 %
## 7.43 %
## 7.53 %
## 7.62 %
## 7.72 %
## 7.82 %
## 7.92 %
## 8.02 %
## 8.11 %
## 8.21 %
## 8.31 %
## 8.41 %
## 8.5 %
## 8.6 %
## 8.7 %
## 8.8 %
## 8.9 %
## 8.99 %
## 9.09 %
## 9.19 %
## 9.29 %
## 9.38 %
## 9.48 %
## 9.58 %
## 9.68 %
## 9.78 %
## 9.87 %
## 9.97 %
## 10.07 %
## 10.17 %
## 10.26 %
## 10.36 %
## 10.46 %
## 10.56 %
## 10.65 %
## 10.75 %
## 10.85 %
## 10.95 %
## 11.05 %
## 11.14 %
## 11.24 %
## 11.34 %
## 11.44 %
## 11.53 %
## 11.63 %
## 11.73 %
## 11.83 %
## 11.93 %
## 12.02 %
## 12.12 %
## 12.22 %
## 12.32 %
## 12.41 %
## 12.51 %
## 12.61 %
## 12.71 %
## 12.81 %
## 12.9 %
## 13 %
## 13.1 %
## 13.2 %
## 13.29 %
## 13.39 %
## 13.49 %
## 13.59 %
## 13.69 %
## 13.78 %
## 13.88 %
## 13.98 %
## 14.08 %
## 14.17 %
## 14.27 %
## 14.37 %
## 14.47 %
## 14.57 %
## 14.66 %
## 14.76 %
## 14.86 %
## 14.96 %
## 15.05 %
## 15.15 %
## 15.25 %
## 15.35 %
## 15.44 %
## 15.54 %
## 15.64 %
## 15.74 %
## 15.84 %
## 15.93 %
## 16.03 %
## 16.13 %
## 16.23 %
## 16.32 %
## 16.42 %
## 16.52 %
## 16.62 %
## 16.72 %
## 16.81 %
## 16.91 %
## 17.01 %
## 17.11 %
## 17.2 %
## 17.3 %
## 17.4 %
## 17.5 %
## 17.6 %
## 17.69 %
## 17.79 %
## 17.89 %
## 17.99 %
## 18.08 %
## 18.18 %
## 18.28 %
## 18.38 %
## 18.48 %
## 18.57 %
## 18.67 %
## 18.77 %
## 18.87 %
## 18.96 %
## 19.06 %
## 19.16 %
## 19.26 %
## 19.35 %
## 19.45 %
## 19.55 %
## 19.65 %
## 19.75 %
## 19.84 %
## 19.94 %
## 20.04 %
## 20.14 %
## 20.23 %
## 20.33 %
## 20.43 %
## 20.53 %
## 20.63 %
## 20.72 %
## 20.82 %
## 20.92 %
## 21.02 %
## 21.11 %
## 21.21 %
## 21.31 %
## 21.41 %
## 21.51 %
## 21.6 %
## 21.7 %
## 21.8 %
## 21.9 %
## 21.99 %
## 22.09 %
## 22.19 %
## 22.29 %
## 22.39 %
## 22.48 %
## 22.58 %
## 22.68 %
## 22.78 %
## 22.87 %
## 22.97 %
## 23.07 %
## 23.17 %
## 23.26 %
## 23.36 %
## 23.46 %
## 23.56 %
## 23.66 %
## 23.75 %
## 23.85 %
## 23.95 %
## 24.05 %
## 24.14 %
## 24.24 %
## 24.34 %
## 24.44 %
## 24.54 %
## 24.63 %
## 24.73 %
## 24.83 %
## 24.93 %
## 25.02 %
## 25.12 %
## 25.22 %
## 25.32 %
## 25.42 %
## 25.51 %
## 25.61 %
## 25.71 %
## 25.81 %
## 25.9 %
## 26 %
## 26.1 %
## 26.2 %
## 26.3 %
## 26.39 %
## 26.49 %
## 26.59 %
## 26.69 %
## 26.78 %
## 26.88 %
## 26.98 %
## 27.08 %
## 27.17 %
## 27.27 %
## 27.37 %
## 27.47 %
## 27.57 %
## 27.66 %
## 27.76 %
## 27.86 %
## 27.96 %
## 28.05 %
## 28.15 %
## 28.25 %
## 28.35 %
## 28.45 %
## 28.54 %
## 28.64 %
## 28.74 %
## 28.84 %
## 28.93 %
## 29.03 %
## 29.13 %
## 29.23 %
## 29.33 %
## 29.42 %
## 29.52 %
## 29.62 %
## 29.72 %
## 29.81 %
## 29.91 %
## 30.01 %
## 30.11 %
## 30.21 %
## 30.3 %
## 30.4 %
## 30.5 %
## 30.6 %
## 30.69 %
## 30.79 %
## 30.89 %
## 30.99 %
## 31.09 %
## 31.18 %
## 31.28 %
## 31.38 %
## 31.48 %
## 31.57 %
## 31.67 %
## 31.77 %
## 31.87 %
## 31.96 %
## 32.06 %
## 32.16 %
## 32.26 %
## 32.36 %
## 32.45 %
## 32.55 %
## 32.65 %
## 32.75 %
## 32.84 %
## 32.94 %
## 33.04 %
## 33.14 %
## 33.24 %
## 33.33 %
## 33.43 %
## 33.53 %
## 33.63 %
## 33.72 %
## 33.82 %
## 33.92 %
## 34.02 %
## 34.12 %
## 34.21 %
## 34.31 %
## 34.41 %
## 34.51 %
## 34.6 %
## 34.7 %
## 34.8 %
## 34.9 %
## 35 %
## 35.09 %
## 35.19 %
## 35.29 %
## 35.39 %
## 35.48 %
## 35.58 %
## 35.68 %
## 35.78 %
## 35.87 %
## 35.97 %
## 36.07 %
## 36.17 %
## 36.27 %
## 36.36 %
## 36.46 %
## 36.56 %
## 36.66 %
## 36.75 %
## 36.85 %
## 36.95 %
## 37.05 %
## 37.15 %
## 37.24 %
## 37.34 %
## 37.44 %
## 37.54 %
## 37.63 %
## 37.73 %
## 37.83 %
## 37.93 %
## 38.03 %
## 38.12 %
## 38.22 %
## 38.32 %
## 38.42 %
## 38.51 %
## 38.61 %
## 38.71 %
## 38.81 %
## 38.91 %
## 39 %
## 39.1 %
## 39.2 %
## 39.3 %
## 39.39 %
## 39.49 %
## 39.59 %
## 39.69 %
## 39.78 %
## 39.88 %
## 39.98 %
## 40.08 %
## 40.18 %
## 40.27 %
## 40.37 %
## 40.47 %
## 40.57 %
## 40.66 %
## 40.76 %
## 40.86 %
## 40.96 %
## 41.06 %
## 41.15 %
## 41.25 %
## 41.35 %
## 41.45 %
## 41.54 %
## 41.64 %
## 41.74 %
## 41.84 %
## 41.94 %
## 42.03 %
## 42.13 %
## 42.23 %
## 42.33 %
## 42.42 %
## 42.52 %
## 42.62 %
## 42.72 %
## 42.82 %
## 42.91 %
## 43.01 %
## 43.11 %
## 43.21 %
## 43.3 %
## 43.4 %
## 43.5 %
## 43.6 %
## 43.7 %
## 43.79 %
## 43.89 %
## 43.99 %
## 44.09 %
## 44.18 %
## 44.28 %
## 44.38 %
## 44.48 %
## 44.57 %
## 44.67 %
## 44.77 %
## 44.87 %
## 44.97 %
## 45.06 %
## 45.16 %
## 45.26 %
## 45.36 %
## 45.45 %
## 45.55 %
## 45.65 %
## 45.75 %
## 45.85 %
## 45.94 %
## 46.04 %
## 46.14 %
## 46.24 %
## 46.33 %
## 46.43 %
## 46.53 %
## 46.63 %
## 46.73 %
## 46.82 %
## 46.92 %
## 47.02 %
## 47.12 %
## 47.21 %
## 47.31 %
## 47.41 %
## 47.51 %
## 47.61 %
## 47.7 %
## 47.8 %
## 47.9 %
## 48 %
## 48.09 %
## 48.19 %
## 48.29 %
## 48.39 %
## 48.48 %
## 48.58 %
## 48.68 %
## 48.78 %
## 48.88 %
## 48.97 %
## 49.07 %
## 49.17 %
## 49.27 %
## 49.36 %
## 49.46 %
## 49.56 %
## 49.66 %
## 49.76 %
## 49.85 %
## 49.95 %
## 50.05 %
## 50.15 %
## 50.24 %
## 50.34 %
## 50.44 %
## 50.54 %
## 50.64 %
## 50.73 %
## 50.83 %
## 50.93 %
## 51.03 %
## 51.12 %
## 51.22 %
## 51.32 %
## 51.42 %
## 51.52 %
## 51.61 %
## 51.71 %
## 51.81 %
## 51.91 %
## 52 %
## 52.1 %
## 52.2 %
## 52.3 %
## 52.39 %
## 52.49 %
## 52.59 %
## 52.69 %
## 52.79 %
## 52.88 %
## 52.98 %
## 53.08 %
## 53.18 %
## 53.27 %
## 53.37 %
## 53.47 %
## 53.57 %
## 53.67 %
## 53.76 %
## 53.86 %
## 53.96 %
## 54.06 %
## 54.15 %
## 54.25 %
## 54.35 %
## 54.45 %
## 54.55 %
## 54.64 %
## 54.74 %
## 54.84 %
## 54.94 %
## 55.03 %
## 55.13 %
## 55.23 %
## 55.33 %
## 55.43 %
## 55.52 %
## 55.62 %
## 55.72 %
## 55.82 %
## 55.91 %
## 56.01 %
## 56.11 %
## 56.21 %
## 56.3 %
## 56.4 %
## 56.5 %
## 56.6 %
## 56.7 %
## 56.79 %
## 56.89 %
## 56.99 %
## 57.09 %
## 57.18 %
## 57.28 %
## 57.38 %
## 57.48 %
## 57.58 %
## 57.67 %
## 57.77 %
## 57.87 %
## 57.97 %
## 58.06 %
## 58.16 %
## 58.26 %
## 58.36 %
## 58.46 %
## 58.55 %
## 58.65 %
## 58.75 %
## 58.85 %
## 58.94 %
## 59.04 %
## 59.14 %
## 59.24 %
## 59.34 %
## 59.43 %
## 59.53 %
## 59.63 %
## 59.73 %
## 59.82 %
## 59.92 %
## 60.02 %
## 60.12 %
## 60.22 %
## 60.31 %
## 60.41 %
## 60.51 %
## 60.61 %
## 60.7 %
## 60.8 %
## 60.9 %
## 61 %
## 61.09 %
## 61.19 %
## 61.29 %
## 61.39 %
## 61.49 %
## 61.58 %
## 61.68 %
## 61.78 %
## 61.88 %
## 61.97 %
## 62.07 %
## 62.17 %
## 62.27 %
## 62.37 %
## 62.46 %
## 62.56 %
## 62.66 %
## 62.76 %
## 62.85 %
## 62.95 %
## 63.05 %
## 63.15 %
## 63.25 %
## 63.34 %
## 63.44 %
## 63.54 %
## 63.64 %
## 63.73 %
## 63.83 %
## 63.93 %
## 64.03 %
## 64.13 %
## 64.22 %
## 64.32 %
## 64.42 %
## 64.52 %
## 64.61 %
## 64.71 %
## 64.81 %
## 64.91 %
## 65 %
## 65.1 %
## 65.2 %
## 65.3 %
## 65.4 %
## 65.49 %
## 65.59 %
## 65.69 %
## 65.79 %
## 65.88 %
## 65.98 %
## 66.08 %
## 66.18 %
## 66.28 %
## 66.37 %
## 66.47 %
## 66.57 %
## 66.67 %
## 66.76 %
## 66.86 %
## 66.96 %
## 67.06 %
## 67.16 %
## 67.25 %
## 67.35 %
## 67.45 %
## 67.55 %
## 67.64 %
## 67.74 %
## 67.84 %
## 67.94 %
## 68.04 %
## 68.13 %
## 68.23 %
## 68.33 %
## 68.43 %
## 68.52 %
## 68.62 %
## 68.72 %
## 68.82 %
## 68.91 %
## 69.01 %
## 69.11 %
## 69.21 %
## 69.31 %
## 69.4 %
## 69.5 %
## 69.6 %
## 69.7 %
## 69.79 %
## 69.89 %
## 69.99 %
## 70.09 %
## 70.19 %
## 70.28 %
## 70.38 %
## 70.48 %
## 70.58 %
## 70.67 %
## 70.77 %
## 70.87 %
## 70.97 %
## 71.07 %
## 71.16 %
## 71.26 %
## 71.36 %
## 71.46 %
## 71.55 %
## 71.65 %
## 71.75 %
## 71.85 %
## 71.95 %
## 72.04 %
## 72.14 %
## 72.24 %
## 72.34 %
## 72.43 %
## 72.53 %
## 72.63 %
## 72.73 %
## 72.83 %
## 72.92 %
## 73.02 %
## 73.12 %
## 73.22 %
## 73.31 %
## 73.41 %
## 73.51 %
## 73.61 %
## 73.7 %
## 73.8 %
## 73.9 %
## 74 %
## 74.1 %
## 74.19 %
## 74.29 %
## 74.39 %
## 74.49 %
## 74.58 %
## 74.68 %
## 74.78 %
## 74.88 %
## 74.98 %
## 75.07 %
## 75.17 %
## 75.27 %
## 75.37 %
## 75.46 %
## 75.56 %
## 75.66 %
## 75.76 %
## 75.86 %
## 75.95 %
## 76.05 %
## 76.15 %
## 76.25 %
## 76.34 %
## 76.44 %
## 76.54 %
## 76.64 %
## 76.74 %
## 76.83 %
## 76.93 %
## 77.03 %
## 77.13 %
## 77.22 %
## 77.32 %
## 77.42 %
## 77.52 %
## 77.61 %
## 77.71 %
## 77.81 %
## 77.91 %
## 78.01 %
## 78.1 %
## 78.2 %
## 78.3 %
## 78.4 %
## 78.49 %
## 78.59 %
## 78.69 %
## 78.79 %
## 78.89 %
## 78.98 %
## 79.08 %
## 79.18 %
## 79.28 %
## 79.37 %
## 79.47 %
## 79.57 %
## 79.67 %
## 79.77 %
## 79.86 %
## 79.96 %
## 80.06 %
## 80.16 %
## 80.25 %
## 80.35 %
## 80.45 %
## 80.55 %
## 80.65 %
## 80.74 %
## 80.84 %
## 80.94 %
## 81.04 %
## 81.13 %
## 81.23 %
## 81.33 %
## 81.43 %
## 81.52 %
## 81.62 %
## 81.72 %
## 81.82 %
## 81.92 %
## 82.01 %
## 82.11 %
## 82.21 %
## 82.31 %
## 82.4 %
## 82.5 %
## 82.6 %
## 82.7 %
## 82.8 %
## 82.89 %
## 82.99 %
## 83.09 %
## 83.19 %
## 83.28 %
## 83.38 %
## 83.48 %
## 83.58 %
## 83.68 %
## 83.77 %
## 83.87 %
## 83.97 %
## 84.07 %
## 84.16 %
## 84.26 %
## 84.36 %
## 84.46 %
## 84.56 %
## 84.65 %
## 84.75 %
## 84.85 %
## 84.95 %
## 85.04 %
## 85.14 %
## 85.24 %
## 85.34 %
## 85.43 %
## 85.53 %
## 85.63 %
## 85.73 %
## 85.83 %
## 85.92 %
## 86.02 %
## 86.12 %
## 86.22 %
## 86.31 %
## 86.41 %
## 86.51 %
## 86.61 %
## 86.71 %
## 86.8 %
## 86.9 %
## 87 %
## 87.1 %
## 87.19 %
## 87.29 %
## 87.39 %
## 87.49 %
## 87.59 %
## 87.68 %
## 87.78 %
## 87.88 %
## 87.98 %
## 88.07 %
## 88.17 %
## 88.27 %
## 88.37 %
## 88.47 %
## 88.56 %
## 88.66 %
## 88.76 %
## 88.86 %
## 88.95 %
## 89.05 %
## 89.15 %
## 89.25 %
## 89.35 %
## 89.44 %
## 89.54 %
## 89.64 %
## 89.74 %
## 89.83 %
## 89.93 %
## 90.03 %
## 90.13 %
## 90.22 %
## 90.32 %
## 90.42 %
## 90.52 %
## 90.62 %
## 90.71 %
## 90.81 %
## 90.91 %
## 91.01 %
## 91.1 %
## 91.2 %
## 91.3 %
## 91.4 %
## 91.5 %
## 91.59 %
## 91.69 %
## 91.79 %
## 91.89 %
## 91.98 %
## 92.08 %
## 92.18 %
## 92.28 %
## 92.38 %
## 92.47 %
## 92.57 %
## 92.67 %
## 92.77 %
## 92.86 %
## 92.96 %
## 93.06 %
## 93.16 %
## 93.26 %
## 93.35 %
## 93.45 %
## 93.55 %
## 93.65 %
## 93.74 %
## 93.84 %
## 93.94 %
## 94.04 %
## 94.13 %
## 94.23 %
## 94.33 %
## 94.43 %
## 94.53 %
## 94.62 %
## 94.72 %
## 94.82 %
## 94.92 %
## 95.01 %
## 95.11 %
## 95.21 %
## 95.31 %
## 95.41 %
## 95.5 %
## 95.6 %
## 95.7 %
## 95.8 %
## 95.89 %
## 95.99 %
## 96.09 %
## 96.19 %
## 96.29 %
## 96.38 %
## 96.48 %
## 96.58 %
## 96.68 %
## 96.77 %
## 96.87 %
## 96.97 %
## 97.07 %
## 97.17 %
## 97.26 %
## 97.36 %
## 97.46 %
## 97.56 %
## 97.65 %
## 97.75 %
## 97.85 %
## 97.95 %
## 98.04 %
## 98.14 %
## 98.24 %
## 98.34 %
## 98.44 %
## 98.53 %
## 98.63 %
## 98.73 %
## 98.83 %
## 98.92 %
## 99.02 %
## 99.12 %
## 99.22 %
## 99.32 %
## 99.41 %
## 99.51 %
## 99.61 %
## 99.71 %
## 99.8 %
## 99.9 %
## 100 %
## Done.
## Creating graphs... Done.
output$plot.p
output$plot.est
# Run glm.sim with graphs, best model, and quadratic terms
output <- glm.sim(data = africa.subset, dv = "gdp_pc", graphs = T, graph.best = 1, 
    quadratics = T)
## Creating model set...
## Number of model terms: 1... 2... 3... 4... 5... 6... 7... 8... Done.
## Running glm...
## Percent complete: 0.39 %
## 0.78 %
## 1.18 %
## 1.57 %
## 1.96 %
## 2.35 %
## 2.75 %
## 3.14 %
## 3.53 %
## 3.92 %
## 4.31 %
## 4.71 %
## 5.1 %
## 5.49 %
## 5.88 %
## 6.27 %
## 6.67 %
## 7.06 %
## 7.45 %
## 7.84 %
## 8.24 %
## 8.63 %
## 9.02 %
## 9.41 %
## 9.8 %
## 10.2 %
## 10.59 %
## 10.98 %
## 11.37 %
## 11.76 %
## 12.16 %
## 12.55 %
## 12.94 %
## 13.33 %
## 13.73 %
## 14.12 %
## 14.51 %
## 14.9 %
## 15.29 %
## 15.69 %
## 16.08 %
## 16.47 %
## 16.86 %
## 17.25 %
## 17.65 %
## 18.04 %
## 18.43 %
## 18.82 %
## 19.22 %
## 19.61 %
## 20 %
## 20.39 %
## 20.78 %
## 21.18 %
## 21.57 %
## 21.96 %
## 22.35 %
## 22.75 %
## 23.14 %
## 23.53 %
## 23.92 %
## 24.31 %
## 24.71 %
## 25.1 %
## 25.49 %
## 25.88 %
## 26.27 %
## 26.67 %
## 27.06 %
## 27.45 %
## 27.84 %
## 28.24 %
## 28.63 %
## 29.02 %
## 29.41 %
## 29.8 %
## 30.2 %
## 30.59 %
## 30.98 %
## 31.37 %
## 31.76 %
## 32.16 %
## 32.55 %
## 32.94 %
## 33.33 %
## 33.73 %
## 34.12 %
## 34.51 %
## 34.9 %
## 35.29 %
## 35.69 %
## 36.08 %
## 36.47 %
## 36.86 %
## 37.25 %
## 37.65 %
## 38.04 %
## 38.43 %
## 38.82 %
## 39.22 %
## 39.61 %
## 40 %
## 40.39 %
## 40.78 %
## 41.18 %
## 41.57 %
## 41.96 %
## 42.35 %
## 42.75 %
## 43.14 %
## 43.53 %
## 43.92 %
## 44.31 %
## 44.71 %
## 45.1 %
## 45.49 %
## 45.88 %
## 46.27 %
## 46.67 %
## 47.06 %
## 47.45 %
## 47.84 %
## 48.24 %
## 48.63 %
## 49.02 %
## 49.41 %
## 49.8 %
## 50.2 %
## 50.59 %
## 50.98 %
## 51.37 %
## 51.76 %
## 52.16 %
## 52.55 %
## 52.94 %
## 53.33 %
## 53.73 %
## 54.12 %
## 54.51 %
## 54.9 %
## 55.29 %
## 55.69 %
## 56.08 %
## 56.47 %
## 56.86 %
## 57.25 %
## 57.65 %
## 58.04 %
## 58.43 %
## 58.82 %
## 59.22 %
## 59.61 %
## 60 %
## 60.39 %
## 60.78 %
## 61.18 %
## 61.57 %
## 61.96 %
## 62.35 %
## 62.75 %
## 63.14 %
## 63.53 %
## 63.92 %
## 64.31 %
## 64.71 %
## 65.1 %
## 65.49 %
## 65.88 %
## 66.27 %
## 66.67 %
## 67.06 %
## 67.45 %
## 67.84 %
## 68.24 %
## 68.63 %
## 69.02 %
## 69.41 %
## 69.8 %
## 70.2 %
## 70.59 %
## 70.98 %
## 71.37 %
## 71.76 %
## 72.16 %
## 72.55 %
## 72.94 %
## 73.33 %
## 73.73 %
## 74.12 %
## 74.51 %
## 74.9 %
## 75.29 %
## 75.69 %
## 76.08 %
## 76.47 %
## 76.86 %
## 77.25 %
## 77.65 %
## 78.04 %
## 78.43 %
## 78.82 %
## 79.22 %
## 79.61 %
## 80 %
## 80.39 %
## 80.78 %
## 81.18 %
## 81.57 %
## 81.96 %
## 82.35 %
## 82.75 %
## 83.14 %
## 83.53 %
## 83.92 %
## 84.31 %
## 84.71 %
## 85.1 %
## 85.49 %
## 85.88 %
## 86.27 %
## 86.67 %
## 87.06 %
## 87.45 %
## 87.84 %
## 88.24 %
## 88.63 %
## 89.02 %
## 89.41 %
## 89.8 %
## 90.2 %
## 90.59 %
## 90.98 %
## 91.37 %
## 91.76 %
## 92.16 %
## 92.55 %
## 92.94 %
## 93.33 %
## 93.73 %
## 94.12 %
## 94.51 %
## 94.9 %
## 95.29 %
## 95.69 %
## 96.08 %
## 96.47 %
## 96.86 %
## 97.25 %
## 97.65 %
## 98.04 %
## 98.43 %
## 98.82 %
## 99.22 %
## 99.61 %
## 100 %
## Done.
## Identifying best models (est)...
## 0.39 %
## 0.78 %
## 1.18 %
## 1.57 %
## 1.96 %
## 2.35 %
## 2.75 %
## 3.14 %
## 3.53 %
## 3.92 %
## 4.31 %
## 4.71 %
## 5.1 %
## 5.49 %
## 5.88 %
## 6.27 %
## 6.67 %
## 7.06 %
## 7.45 %
## 7.84 %
## 8.24 %
## 8.63 %
## 9.02 %
## 9.41 %
## 9.8 %
## 10.2 %
## 10.59 %
## 10.98 %
## 11.37 %
## 11.76 %
## 12.16 %
## 12.55 %
## 12.94 %
## 13.33 %
## 13.73 %
## 14.12 %
## 14.51 %
## 14.9 %
## 15.29 %
## 15.69 %
## 16.08 %
## 16.47 %
## 16.86 %
## 17.25 %
## 17.65 %
## 18.04 %
## 18.43 %
## 18.82 %
## 19.22 %
## 19.61 %
## 20 %
## 20.39 %
## 20.78 %
## 21.18 %
## 21.57 %
## 21.96 %
## 22.35 %
## 22.75 %
## 23.14 %
## 23.53 %
## 23.92 %
## 24.31 %
## 24.71 %
## 25.1 %
## 25.49 %
## 25.88 %
## 26.27 %
## 26.67 %
## 27.06 %
## 27.45 %
## 27.84 %
## 28.24 %
## 28.63 %
## 29.02 %
## 29.41 %
## 29.8 %
## 30.2 %
## 30.59 %
## 30.98 %
## 31.37 %
## 31.76 %
## 32.16 %
## 32.55 %
## 32.94 %
## 33.33 %
## 33.73 %
## 34.12 %
## 34.51 %
## 34.9 %
## 35.29 %
## 35.69 %
## 36.08 %
## 36.47 %
## 36.86 %
## 37.25 %
## 37.65 %
## 38.04 %
## 38.43 %
## 38.82 %
## 39.22 %
## 39.61 %
## 40 %
## 40.39 %
## 40.78 %
## 41.18 %
## 41.57 %
## 41.96 %
## 42.35 %
## 42.75 %
## 43.14 %
## 43.53 %
## 43.92 %
## 44.31 %
## 44.71 %
## 45.1 %
## 45.49 %
## 45.88 %
## 46.27 %
## 46.67 %
## 47.06 %
## 47.45 %
## 47.84 %
## 48.24 %
## 48.63 %
## 49.02 %
## 49.41 %
## 49.8 %
## 50.2 %
## 50.59 %
## 50.98 %
## 51.37 %
## 51.76 %
## 52.16 %
## 52.55 %
## 52.94 %
## 53.33 %
## 53.73 %
## 54.12 %
## 54.51 %
## 54.9 %
## 55.29 %
## 55.69 %
## 56.08 %
## 56.47 %
## 56.86 %
## 57.25 %
## 57.65 %
## 58.04 %
## 58.43 %
## 58.82 %
## 59.22 %
## 59.61 %
## 60 %
## 60.39 %
## 60.78 %
## 61.18 %
## 61.57 %
## 61.96 %
## 62.35 %
## 62.75 %
## 63.14 %
## 63.53 %
## 63.92 %
## 64.31 %
## 64.71 %
## 65.1 %
## 65.49 %
## 65.88 %
## 66.27 %
## 66.67 %
## 67.06 %
## 67.45 %
## 67.84 %
## 68.24 %
## 68.63 %
## 69.02 %
## 69.41 %
## 69.8 %
## 70.2 %
## 70.59 %
## 70.98 %
## 71.37 %
## 71.76 %
## 72.16 %
## 72.55 %
## 72.94 %
## 73.33 %
## 73.73 %
## 74.12 %
## 74.51 %
## 74.9 %
## 75.29 %
## 75.69 %
## 76.08 %
## 76.47 %
## 76.86 %
## 77.25 %
## 77.65 %
## 78.04 %
## 78.43 %
## 78.82 %
## 79.22 %
## 79.61 %
## 80 %
## 80.39 %
## 80.78 %
## 81.18 %
## 81.57 %
## 81.96 %
## 82.35 %
## 82.75 %
## 83.14 %
## 83.53 %
## 83.92 %
## 84.31 %
## 84.71 %
## 85.1 %
## 85.49 %
## 85.88 %
## 86.27 %
## 86.67 %
## 87.06 %
## 87.45 %
## 87.84 %
## 88.24 %
## 88.63 %
## 89.02 %
## 89.41 %
## 89.8 %
## 90.2 %
## 90.59 %
## 90.98 %
## 91.37 %
## 91.76 %
## 92.16 %
## 92.55 %
## 92.94 %
## 93.33 %
## 93.73 %
## 94.12 %
## 94.51 %
## 94.9 %
## 95.29 %
## 95.69 %
## 96.08 %
## 96.47 %
## 96.86 %
## 97.25 %
## 97.65 %
## 98.04 %
## 98.43 %
## 98.82 %
## 99.22 %
## 99.61 %
## 100 %
## Binding data...
## Arranging data...
## Identifying best models (p)...
## 0.39 %
## 0.78 %
## 1.18 %
## 1.57 %
## 1.96 %
## 2.35 %
## 2.75 %
## 3.14 %
## 3.53 %
## 3.92 %
## 4.31 %
## 4.71 %
## 5.1 %
## 5.49 %
## 5.88 %
## 6.27 %
## 6.67 %
## 7.06 %
## 7.45 %
## 7.84 %
## 8.24 %
## 8.63 %
## 9.02 %
## 9.41 %
## 9.8 %
## 10.2 %
## 10.59 %
## 10.98 %
## 11.37 %
## 11.76 %
## 12.16 %
## 12.55 %
## 12.94 %
## 13.33 %
## 13.73 %
## 14.12 %
## 14.51 %
## 14.9 %
## 15.29 %
## 15.69 %
## 16.08 %
## 16.47 %
## 16.86 %
## 17.25 %
## 17.65 %
## 18.04 %
## 18.43 %
## 18.82 %
## 19.22 %
## 19.61 %
## 20 %
## 20.39 %
## 20.78 %
## 21.18 %
## 21.57 %
## 21.96 %
## 22.35 %
## 22.75 %
## 23.14 %
## 23.53 %
## 23.92 %
## 24.31 %
## 24.71 %
## 25.1 %
## 25.49 %
## 25.88 %
## 26.27 %
## 26.67 %
## 27.06 %
## 27.45 %
## 27.84 %
## 28.24 %
## 28.63 %
## 29.02 %
## 29.41 %
## 29.8 %
## 30.2 %
## 30.59 %
## 30.98 %
## 31.37 %
## 31.76 %
## 32.16 %
## 32.55 %
## 32.94 %
## 33.33 %
## 33.73 %
## 34.12 %
## 34.51 %
## 34.9 %
## 35.29 %
## 35.69 %
## 36.08 %
## 36.47 %
## 36.86 %
## 37.25 %
## 37.65 %
## 38.04 %
## 38.43 %
## 38.82 %
## 39.22 %
## 39.61 %
## 40 %
## 40.39 %
## 40.78 %
## 41.18 %
## 41.57 %
## 41.96 %
## 42.35 %
## 42.75 %
## 43.14 %
## 43.53 %
## 43.92 %
## 44.31 %
## 44.71 %
## 45.1 %
## 45.49 %
## 45.88 %
## 46.27 %
## 46.67 %
## 47.06 %
## 47.45 %
## 47.84 %
## 48.24 %
## 48.63 %
## 49.02 %
## 49.41 %
## 49.8 %
## 50.2 %
## 50.59 %
## 50.98 %
## 51.37 %
## 51.76 %
## 52.16 %
## 52.55 %
## 52.94 %
## 53.33 %
## 53.73 %
## 54.12 %
## 54.51 %
## 54.9 %
## 55.29 %
## 55.69 %
## 56.08 %
## 56.47 %
## 56.86 %
## 57.25 %
## 57.65 %
## 58.04 %
## 58.43 %
## 58.82 %
## 59.22 %
## 59.61 %
## 60 %
## 60.39 %
## 60.78 %
## 61.18 %
## 61.57 %
## 61.96 %
## 62.35 %
## 62.75 %
## 63.14 %
## 63.53 %
## 63.92 %
## 64.31 %
## 64.71 %
## 65.1 %
## 65.49 %
## 65.88 %
## 66.27 %
## 66.67 %
## 67.06 %
## 67.45 %
## 67.84 %
## 68.24 %
## 68.63 %
## 69.02 %
## 69.41 %
## 69.8 %
## 70.2 %
## 70.59 %
## 70.98 %
## 71.37 %
## 71.76 %
## 72.16 %
## 72.55 %
## 72.94 %
## 73.33 %
## 73.73 %
## 74.12 %
## 74.51 %
## 74.9 %
## 75.29 %
## 75.69 %
## 76.08 %
## 76.47 %
## 76.86 %
## 77.25 %
## 77.65 %
## 78.04 %
## 78.43 %
## 78.82 %
## 79.22 %
## 79.61 %
## 80 %
## 80.39 %
## 80.78 %
## 81.18 %
## 81.57 %
## 81.96 %
## 82.35 %
## 82.75 %
## 83.14 %
## 83.53 %
## 83.92 %
## 84.31 %
## 84.71 %
## 85.1 %
## 85.49 %
## 85.88 %
## 86.27 %
## 86.67 %
## 87.06 %
## 87.45 %
## 87.84 %
## 88.24 %
## 88.63 %
## 89.02 %
## 89.41 %
## 89.8 %
## 90.2 %
## 90.59 %
## 90.98 %
## 91.37 %
## 91.76 %
## 92.16 %
## 92.55 %
## 92.94 %
## 93.33 %
## 93.73 %
## 94.12 %
## 94.51 %
## 94.9 %
## 95.29 %
## 95.69 %
## 96.08 %
## 96.47 %
## 96.86 %
## 97.25 %
## 97.65 %
## 98.04 %
## 98.43 %
## 98.82 %
## 99.22 %
## 99.61 %
## 100 %
## Binding data...
## Arranging data...
## Done.
## Extracting results...
## 0.39 %
## 0.78 %
## 1.18 %
## 1.57 %
## 1.96 %
## 2.35 %
## 2.75 %
## 3.14 %
## 3.53 %
## 3.92 %
## 4.31 %
## 4.71 %
## 5.1 %
## 5.49 %
## 5.88 %
## 6.27 %
## 6.67 %
## 7.06 %
## 7.45 %
## 7.84 %
## 8.24 %
## 8.63 %
## 9.02 %
## 9.41 %
## 9.8 %
## 10.2 %
## 10.59 %
## 10.98 %
## 11.37 %
## 11.76 %
## 12.16 %
## 12.55 %
## 12.94 %
## 13.33 %
## 13.73 %
## 14.12 %
## 14.51 %
## 14.9 %
## 15.29 %
## 15.69 %
## 16.08 %
## 16.47 %
## 16.86 %
## 17.25 %
## 17.65 %
## 18.04 %
## 18.43 %
## 18.82 %
## 19.22 %
## 19.61 %
## 20 %
## 20.39 %
## 20.78 %
## 21.18 %
## 21.57 %
## 21.96 %
## 22.35 %
## 22.75 %
## 23.14 %
## 23.53 %
## 23.92 %
## 24.31 %
## 24.71 %
## 25.1 %
## 25.49 %
## 25.88 %
## 26.27 %
## 26.67 %
## 27.06 %
## 27.45 %
## 27.84 %
## 28.24 %
## 28.63 %
## 29.02 %
## 29.41 %
## 29.8 %
## 30.2 %
## 30.59 %
## 30.98 %
## 31.37 %
## 31.76 %
## 32.16 %
## 32.55 %
## 32.94 %
## 33.33 %
## 33.73 %
## 34.12 %
## 34.51 %
## 34.9 %
## 35.29 %
## 35.69 %
## 36.08 %
## 36.47 %
## 36.86 %
## 37.25 %
## 37.65 %
## 38.04 %
## 38.43 %
## 38.82 %
## 39.22 %
## 39.61 %
## 40 %
## 40.39 %
## 40.78 %
## 41.18 %
## 41.57 %
## 41.96 %
## 42.35 %
## 42.75 %
## 43.14 %
## 43.53 %
## 43.92 %
## 44.31 %
## 44.71 %
## 45.1 %
## 45.49 %
## 45.88 %
## 46.27 %
## 46.67 %
## 47.06 %
## 47.45 %
## 47.84 %
## 48.24 %
## 48.63 %
## 49.02 %
## 49.41 %
## 49.8 %
## 50.2 %
## 50.59 %
## 50.98 %
## 51.37 %
## 51.76 %
## 52.16 %
## 52.55 %
## 52.94 %
## 53.33 %
## 53.73 %
## 54.12 %
## 54.51 %
## 54.9 %
## 55.29 %
## 55.69 %
## 56.08 %
## 56.47 %
## 56.86 %
## 57.25 %
## 57.65 %
## 58.04 %
## 58.43 %
## 58.82 %
## 59.22 %
## 59.61 %
## 60 %
## 60.39 %
## 60.78 %
## 61.18 %
## 61.57 %
## 61.96 %
## 62.35 %
## 62.75 %
## 63.14 %
## 63.53 %
## 63.92 %
## 64.31 %
## 64.71 %
## 65.1 %
## 65.49 %
## 65.88 %
## 66.27 %
## 66.67 %
## 67.06 %
## 67.45 %
## 67.84 %
## 68.24 %
## 68.63 %
## 69.02 %
## 69.41 %
## 69.8 %
## 70.2 %
## 70.59 %
## 70.98 %
## 71.37 %
## 71.76 %
## 72.16 %
## 72.55 %
## 72.94 %
## 73.33 %
## 73.73 %
## 74.12 %
## 74.51 %
## 74.9 %
## 75.29 %
## 75.69 %
## 76.08 %
## 76.47 %
## 76.86 %
## 77.25 %
## 77.65 %
## 78.04 %
## 78.43 %
## 78.82 %
## 79.22 %
## 79.61 %
## 80 %
## 80.39 %
## 80.78 %
## 81.18 %
## 81.57 %
## 81.96 %
## 82.35 %
## 82.75 %
## 83.14 %
## 83.53 %
## 83.92 %
## 84.31 %
## 84.71 %
## 85.1 %
## 85.49 %
## 85.88 %
## 86.27 %
## 86.67 %
## 87.06 %
## 87.45 %
## 87.84 %
## 88.24 %
## 88.63 %
## 89.02 %
## 89.41 %
## 89.8 %
## 90.2 %
## 90.59 %
## 90.98 %
## 91.37 %
## 91.76 %
## 92.16 %
## 92.55 %
## 92.94 %
## 93.33 %
## 93.73 %
## 94.12 %
## 94.51 %
## 94.9 %
## 95.29 %
## 95.69 %
## 96.08 %
## 96.47 %
## 96.86 %
## 97.25 %
## 97.65 %
## 98.04 %
## 98.43 %
## 98.82 %
## 99.22 %
## 99.61 %
## 100 %
## Done.
## Creating graphs... Done.
output$plot.p
output$plot.est

Inferences/Regression Diagnostics

bptest(oldman.reg)
## 
##  studentized Breusch-Pagan test
## 
## data:  oldman.reg
## BP = 1530, df = 14, p-value < 2.2e-16
bgtest(oldman.reg)
## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  oldman.reg
## LM test = 10.35, df = 1, p-value = 0.001292
dwtest(oldman.reg)
## 
##  Durbin-Watson test
## 
## data:  oldman.reg
## DW = 1.975, p-value = 0.0006415
## alternative hypothesis: true autocorrelation is greater than 0
coeftest(x = oldman.reg, vcov. = vcovHC)
## 
## z test of coefficients:
## 
##                        Estimate Std. Error z value Pr(>|z|)    
## (Intercept)             -2.0125     0.0464  -43.42  < 2e-16 ***
## partyidrepublican        0.0951     0.0444    2.14  0.03221 *  
## partyidindependent       0.1283     0.0491    2.61  0.00897 ** 
## partyidsomething else   -0.7358     0.1365   -5.39  7.1e-08 ***
## income$15,000-$29,999   -0.1543     0.0488   -3.16  0.00156 ** 
## income$30,000-$49,999   -0.5768     0.0533  -10.83  < 2e-16 ***
## income$50,000-$74,999   -0.8934     0.0545  -16.39  < 2e-16 ***
## income$75,000-$99,999   -0.8447     0.1190   -7.10  1.3e-12 ***
## relign8catholic         -0.5387     0.0455  -11.83  < 2e-16 ***
## relign8mormon/lds       -0.5416     0.1057   -5.13  3.0e-07 ***
## relign8other christian  -1.8544     0.0726  -25.53  < 2e-16 ***
## relign8jewish            0.1297     0.1058    1.23  0.22033    
## relign8muslim           -1.7141     0.4556   -3.76  0.00017 ***
## relign8something else   -1.7263     0.1232  -14.01  < 2e-16 ***
## relign8none             -1.1830     0.0783  -15.11  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

By way of introduction...

quick/dirty v. involved/fancy

tips <- reshape2::tips  # Load dataset on tipping behavior included with reshape2 package
attributes(tips)  # Check attributes of the tips dataset (names, row.names, class)
## $names
## [1] "total_bill" "tip"        "sex"        "smoker"     "day"       
## [6] "time"       "size"      
## 
## $row.names
##   [1] "1"   "2"   "3"   "4"   "5"   "6"   "7"   "8"   "9"   "10"  "11" 
##  [12] "12"  "13"  "14"  "15"  "16"  "17"  "18"  "19"  "20"  "21"  "22" 
##  [23] "23"  "24"  "25"  "26"  "27"  "28"  "29"  "30"  "31"  "32"  "33" 
##  [34] "34"  "35"  "36"  "37"  "38"  "39"  "40"  "41"  "42"  "43"  "44" 
##  [45] "45"  "46"  "47"  "48"  "49"  "50"  "51"  "52"  "53"  "54"  "55" 
##  [56] "56"  "57"  "58"  "59"  "60"  "61"  "62"  "63"  "64"  "65"  "66" 
##  [67] "67"  "68"  "69"  "70"  "71"  "72"  "73"  "74"  "75"  "76"  "77" 
##  [78] "78"  "79"  "80"  "81"  "82"  "83"  "84"  "85"  "86"  "87"  "88" 
##  [89] "89"  "90"  "91"  "92"  "93"  "94"  "95"  "96"  "97"  "98"  "99" 
## [100] "100" "101" "102" "103" "104" "105" "106" "107" "108" "109" "110"
## [111] "111" "112" "113" "114" "115" "116" "117" "118" "119" "120" "121"
## [122] "122" "123" "124" "125" "126" "127" "128" "129" "130" "131" "132"
## [133] "133" "134" "135" "136" "137" "138" "139" "140" "141" "142" "143"
## [144] "144" "145" "146" "147" "148" "149" "150" "151" "152" "153" "154"
## [155] "155" "156" "157" "158" "159" "160" "161" "162" "163" "164" "165"
## [166] "166" "167" "168" "169" "170" "171" "172" "173" "174" "175" "176"
## [177] "177" "178" "179" "180" "181" "182" "183" "184" "185" "186" "187"
## [188] "188" "189" "190" "191" "192" "193" "194" "195" "196" "197" "198"
## [199] "199" "200" "201" "202" "203" "204" "205" "206" "207" "208" "209"
## [210] "210" "211" "212" "213" "214" "215" "216" "217" "218" "219" "220"
## [221] "221" "222" "223" "224" "225" "226" "227" "228" "229" "230" "231"
## [232] "232" "233" "234" "235" "236" "237" "238" "239" "240" "241" "242"
## [243] "243" "244"
## 
## $class
## [1] "data.frame"
# Create an object of class 'lm' (linear model), regressing tip on some
# covariates
tips.reg <- lm(formula = tip ~ total_bill + sex + smoker + day + time + size, 
    data = tips)
attributes(tips.reg)  # Check attributes of the tips.reg object (names, class)
## $names
##  [1] "coefficients"  "residuals"     "effects"       "rank"         
##  [5] "fitted.values" "assign"        "qr"            "df.residual"  
##  [9] "contrasts"     "xlevels"       "call"          "terms"        
## [13] "model"        
## 
## $class
## [1] "lm"
plot(tips)  # Calls plotting method for class of tips dataset ('data.frame')
plot(tips.reg, which = 1:2)  # Calls plotting method for class of tips.reg objects ('lm'), print first two plots only

xyplot(tips)  # Attempt in lattice to automatically plot objects of class 'data.frame'
## Error: no applicable method for 'xyplot' applied to an object of class
## "data.frame"
ggplot(data = tips) + geom_point()  # Attempt in ggplot to automatically plot objects of class 'data.frame'
## Error: 'where' is missing
xyplot(tips.reg)  # Attempt in lattice to automatically plot objects of class 'lm'
## Error: no applicable method for 'xyplot' applied to an object of class
## "lm"
ggplot(data = tips.reg) + geom_point()  # Attempt in ggplot to automatically plot objects of class 'lm'
## Error: 'where' is missing

The Dataset

For more info: http://www.ipw.unibe.ch/content/team/klaus_armingeon/comparative_political_data_sets/index_eng.html

A copy of the dataset is available in the file "cpds.csv", which can be read in as:

data <- read.csv(file = "datasets/cpds.csv")

The Dataset, more specifically

attributes(data)[1:2]  # Only show the first two attributes of the dataset (column names and object class)
## $names
##   [1] "year"         "country"      "countryn"     "gov_right"   
##   [5] "gov_cent"     "gov_left"     "govparty"     "gov_new"     
##   [9] "gov_gap"      "gov_type"     "gov_chan"     "elect"       
##  [13] "vturn"        "social1"      "social2"      "social3"     
##  [17] "social4"      "leftsoc1"     "leftsoc2"     "leftsoc3"    
##  [21] "comm1"        "comm2"        "comm3"        "comm4"       
##  [25] "agrarian"     "conserv1"     "conserv2"     "conserv3"    
##  [29] "conserv4"     "conserv5"     "relig1"       "relig2"      
##  [33] "relig3"       "relig4"       "relig5"       "relig6"      
##  [37] "relig7"       "liberal1"     "liberal2"     "liberal3"    
##  [41] "liberal4"     "liberal5"     "liberal6"     "ultrar1"     
##  [45] "ultrar2"      "protest1"     "protest2"     "protest3"    
##  [49] "protest4"     "green1"       "green2"       "green3"      
##  [53] "ethnic1"      "ethnic2"      "ethnic3"      "ethnic4"     
##  [57] "leftall"      "centall"      "rightall"     "others"      
##  [61] "VAR00017"     "VAR00019"     "ssocial1"     "ssocial2"    
##  [65] "ssocial3"     "ssocial4"     "sleftso1"     "sleftso2"    
##  [69] "sleftso3"     "scomm1"       "scomm2"       "scomm3"      
##  [73] "scomm4"       "sagraria"     "sconser1"     "sconser2"    
##  [77] "sconser3"     "sconser4"     "sconser5"     "srelig1"     
##  [81] "srelig2"      "srelig3"      "srelig4"      "srelig5"     
##  [85] "srelig6"      "srelig7"      "slibera1"     "slibera2"    
##  [89] "slibera3"     "slibera4"     "slibera5"     "slibera6"    
##  [93] "sultrar1"     "sultrar2"     "sprotes1"     "sprotes2"    
##  [97] "sprotes3"     "sprotes4"     "sgreen1"      "sgreen2"     
## [101] "sgreen3"      "sethnic1"     "sethnic2"     "sethnic3"    
## [105] "sethnic4"     "sleftal"      "scental"      "srightal"    
## [109] "sothers"      "VAR00018"     "VAR00020"     "womenpar"    
## [113] "rae_ele"      "rae_leg"      "effpar_ele"   "effpar_leg"  
## [117] "dis_abso"     "dis_rel"      "dis_gall"     "var00001"    
## [121] "var00002"     "lfirst1"      "lfirst2"      "lsec1"       
## [125] "lsec2"        "leff1"        "leff2"        "lmin1"       
## [129] "lmin2"        "lexe1"        "lexe2"        "ldis1"       
## [133] "ldis2"        "lint1"        "lint2"        "lfed1"       
## [137] "lfed2"        "lbic1"        "lbic2"        "lrid1"       
## [141] "lrid2"        "ljud1"        "ljud2"        "lbank1"      
## [145] "lbank2"       "lfirstp"      "lfirstpi"     "lfirstpb"    
## [149] "var00003"     "var00008"     "instcons"     "plural"      
## [153] "structur"     "integr"       "fed"          "pres"        
## [157] "singmemd"     "strbic"       "referen"      "judrev"      
## [161] "VAR00015"     "VAR00016"     "kaopen"       "openc"       
## [165] "outlays"      "receipts"     "realgdpgr"    "nomgdpgr"    
## [169] "inflation"    "debt"         "deficit"      "adjustdef"   
## [173] "interest"     "VAR00013"     "VAR00014"     "ttl_labf"    
## [177] "civ_labf"     "emp_civ"      "labfopar"     "empratio"    
## [181] "emp_ag"       "emp_ind"      "emp_serv"     "emp_un"      
## [185] "unemp"        "st_unemp"     "var00005"     "var00004"    
## [189] "nld"          "wi"           "wdlost"       "strike"      
## [193] "grossu"       "netu"         "ud"           "adjcov"      
## [197] "VAR00027"     "VAR00028"     "sstran"       "socexp_t_pmp"
## [201] "socexp_c_pmp" "socexp_k_pmp" "oldage_pmp"   "survivor_pmp"
## [205] "incapben_pmp" "health_pmp"   "family_pmp"   "almp_pmp"    
## [209] "unemp_pmp"    "housing_pmp"  "othsocx_pmp"  "var00007"    
## [213] "var00006"     "fallow_pmp"   "mpleave_pmp"  "othfam_c_pmp"
## [217] "daycare_pmp"  "othfam_k_pmp" "VAR00047"     "VAR00048"    
## [221] "servadmi_pmp" "training_pmp" "jobrot_pmp"   "incent_pmp"  
## [225] "disabled_pmp" "jobcrea_pmp"  "startup_pmp"  "compen_pmp"  
## [229] "earretir_pmp" "emprot"       "VAR00009"     "VAR00010"    
## [233] "ilo_tot"      "pop"          "pop15_64"     "pop65"       
## [237] "elderly"     
## 
## $class
## [1] "data.frame"
# There are many variables, so it can be helpful to extract their classes
# (mostly to check for factors) via a quick for-loop

classes <- NULL  # Create a placeholder for the class output

# for-loop
for (i in 1:ncol(data)) {
    classes <- c(classes, class(data[, i]))
}
sort(classes)  # Sort the results alphabetically (appears to be some factors, a few integers, and many numerics)
##   [1] "factor"  "factor"  "factor"  "factor"  "factor"  "factor"  "factor" 
##   [8] "factor"  "factor"  "factor"  "factor"  "factor"  "factor"  "integer"
##  [15] "integer" "integer" "integer" "integer" "integer" "integer" "integer"
##  [22] "integer" "integer" "integer" "integer" "integer" "integer" "integer"
##  [29] "integer" "integer" "integer" "integer" "integer" "integer" "integer"
##  [36] "integer" "integer" "integer" "integer" "integer" "integer" "integer"
##  [43] "integer" "integer" "integer" "integer" "numeric" "numeric" "numeric"
##  [50] "numeric" "numeric" "numeric" "numeric" "numeric" "numeric" "numeric"
##  [57] "numeric" "numeric" "numeric" "numeric" "numeric" "numeric" "numeric"
##  [64] "numeric" "numeric" "numeric" "numeric" "numeric" "numeric" "numeric"
##  [71] "numeric" "numeric" "numeric" "numeric" "numeric" "numeric" "numeric"
##  [78] "numeric" "numeric" "numeric" "numeric" "numeric" "numeric" "numeric"
##  [85] "numeric" "numeric" "numeric" "numeric" "numeric" "numeric" "numeric"
##  [92] "numeric" "numeric" "numeric" "numeric" "numeric" "numeric" "numeric"
##  [99] "numeric" "numeric" "numeric" "numeric" "numeric" "numeric" "numeric"
## [106] "numeric" "numeric" "numeric" "numeric" "numeric" "numeric" "numeric"
## [113] "numeric" "numeric" "numeric" "numeric" "numeric" "numeric" "numeric"
## [120] "numeric" "numeric" "numeric" "numeric" "numeric" "numeric" "numeric"
## [127] "numeric" "numeric" "numeric" "numeric" "numeric" "numeric" "numeric"
## [134] "numeric" "numeric" "numeric" "numeric" "numeric" "numeric" "numeric"
## [141] "numeric" "numeric" "numeric" "numeric" "numeric" "numeric" "numeric"
## [148] "numeric" "numeric" "numeric" "numeric" "numeric" "numeric" "numeric"
## [155] "numeric" "numeric" "numeric" "numeric" "numeric" "numeric" "numeric"
## [162] "numeric" "numeric" "numeric" "numeric" "numeric" "numeric" "numeric"
## [169] "numeric" "numeric" "numeric" "numeric" "numeric" "numeric" "numeric"
## [176] "numeric" "numeric" "numeric" "numeric" "numeric" "numeric" "numeric"
## [183] "numeric" "numeric" "numeric" "numeric" "numeric" "numeric" "numeric"
## [190] "numeric" "numeric" "numeric" "numeric" "numeric" "numeric" "numeric"
## [197] "numeric" "numeric" "numeric" "numeric" "numeric" "numeric" "numeric"
## [204] "numeric" "numeric" "numeric" "numeric" "numeric" "numeric" "numeric"
## [211] "numeric" "numeric" "numeric" "numeric" "numeric" "numeric" "numeric"
## [218] "numeric" "numeric" "numeric" "numeric" "numeric" "numeric" "numeric"
## [225] "numeric" "numeric" "numeric" "numeric" "numeric" "numeric" "numeric"
## [232] "numeric" "numeric" "numeric" "numeric" "numeric" "numeric"

base

plot(x = )
plot(x = vturn)  # Attempt to plot variable 'vturn' without vectorization
## Error: object 'vturn' not found
plot(x = data$vturn)  # Attempt to plot variable 'vturn' with '$' vectorization
plot(x = , y = , type = "")

base/type (scatter, line, both)

plot(x = data$year, y = data$realgdpgr, type = "p")  # Plot Year on x-axis and Real GDP Growth on y-axis
plot(x = data$year, y = data$realgdpgr, type = "l")  # Plot Year on x-axis and Real GDP Growth on y-axis
plot(x = data$year, y = data$realgdpgr, type = "b")  # Plot Year on x-axis and Real GDP Growth on y-axis

base/type (histograms, density plots)

hist(x = data$vturn)  # Plot histogram of voter turnout 
hist(x = data$vturn, breaks = 2)  # Plot histogram of voter turnout, with 2 breaks
hist(x = data$vturn, breaks = 50)  # Plot histogram of voter turnout, with 50 breaks
vturn.density <- density(x = data$vturn, na.rm = T)  # Create a density object (NOTE: be sure to remove missing values)
class(vturn.density)  # Check the class of the object
## [1] "density"
vturn.density  # View the contents of the object
## 
## Call:
##  density.default(x = data$vturn, na.rm = T)
## 
## Data: data$vturn (1134 obs.);    Bandwidth 'bw' = 2.826
## 
##        x               y          
##  Min.   : 26.5   Min.   :0.00001  
##  1st Qu.: 46.0   1st Qu.:0.00248  
##  Median : 65.4   Median :0.00626  
##  Mean   : 65.4   Mean   :0.01285  
##  3rd Qu.: 84.8   3rd Qu.:0.02472  
##  Max.   :104.3   Max.   :0.03609
plot(x = vturn.density)  # Plot the density object
plot(x = density(x = data$vturn, bw = 2, na.rm = T))  # Plot the density object, bandwidth of 2
plot(x = density(x = data$vturn, bw = 0.5, na.rm = T))  # Plot the density object, bandwidth of 0.5
plot(x = density(x = data$vturn, bw = 6, na.rm = T))  # Plot the density object, bandwidth of 6

base/options [labeling]

plot(x = , y = , type = "", xlab = "", ylab = "", main = "")
plot(x = data$year, y = data$realgdpgr, type = "p", xlab = "Year", ylab = "Real GDP Growth", 
    main = "This Graph is Great")  # Add labels for axes and overall plot

base/options [axis + size scaling]

plot(x = , y = , type = "", xlim = , ylim = , cex = )
plot(x = data$year, y = data$realgdpgr, type = "p")  # Create a basic plot
# Limit the years (x-axis) to between 1976 and 1991
plot(x = data$year, y = data$realgdpgr, type = "p", xlim = c(1976, 1991))
# Limit the years (x-axis) to between 1976 and 1991, increase point size to
# 2
plot(x = data$year, y = data$realgdpgr, type = "p", xlim = c(1976, 1991), cex = 2)
# Limit the years (x-axis) to between 1976 and 1991, decrease point size to
# 0.5
plot(x = data$year, y = data$realgdpgr, type = "p", xlim = c(1976, 1991), cex = 0.5)

base/options [graphical parameters]

plot(x = , y = , type = "", col = "", pch = , lty = , lwd = )
colors()  # View all elements of the color vector
##   [1] "white"                "aliceblue"            "antiquewhite"        
##   [4] "antiquewhite1"        "antiquewhite2"        "antiquewhite3"       
##   [7] "antiquewhite4"        "aquamarine"           "aquamarine1"         
##  [10] "aquamarine2"          "aquamarine3"          "aquamarine4"         
##  [13] "azure"                "azure1"               "azure2"              
##  [16] "azure3"               "azure4"               "beige"               
##  [19] "bisque"               "bisque1"              "bisque2"             
##  [22] "bisque3"              "bisque4"              "black"               
##  [25] "blanchedalmond"       "blue"                 "blue1"               
##  [28] "blue2"                "blue3"                "blue4"               
##  [31] "blueviolet"           "brown"                "brown1"              
##  [34] "brown2"               "brown3"               "brown4"              
##  [37] "burlywood"            "burlywood1"           "burlywood2"          
##  [40] "burlywood3"           "burlywood4"           "cadetblue"           
##  [43] "cadetblue1"           "cadetblue2"           "cadetblue3"          
##  [46] "cadetblue4"           "chartreuse"           "chartreuse1"         
##  [49] "chartreuse2"          "chartreuse3"          "chartreuse4"         
##  [52] "chocolate"            "chocolate1"           "chocolate2"          
##  [55] "chocolate3"           "chocolate4"           "coral"               
##  [58] "coral1"               "coral2"               "coral3"              
##  [61] "coral4"               "cornflowerblue"       "cornsilk"            
##  [64] "cornsilk1"            "cornsilk2"            "cornsilk3"           
##  [67] "cornsilk4"            "cyan"                 "cyan1"               
##  [70] "cyan2"                "cyan3"                "cyan4"               
##  [73] "darkblue"             "darkcyan"             "darkgoldenrod"       
##  [76] "darkgoldenrod1"       "darkgoldenrod2"       "darkgoldenrod3"      
##  [79] "darkgoldenrod4"       "darkgray"             "darkgreen"           
##  [82] "darkgrey"             "darkkhaki"            "darkmagenta"         
##  [85] "darkolivegreen"       "darkolivegreen1"      "darkolivegreen2"     
##  [88] "darkolivegreen3"      "darkolivegreen4"      "darkorange"          
##  [91] "darkorange1"          "darkorange2"          "darkorange3"         
##  [94] "darkorange4"          "darkorchid"           "darkorchid1"         
##  [97] "darkorchid2"          "darkorchid3"          "darkorchid4"         
## [100] "darkred"              "darksalmon"           "darkseagreen"        
## [103] "darkseagreen1"        "darkseagreen2"        "darkseagreen3"       
## [106] "darkseagreen4"        "darkslateblue"        "darkslategray"       
## [109] "darkslategray1"       "darkslategray2"       "darkslategray3"      
## [112] "darkslategray4"       "darkslategrey"        "darkturquoise"       
## [115] "darkviolet"           "deeppink"             "deeppink1"           
## [118] "deeppink2"            "deeppink3"            "deeppink4"           
## [121] "deepskyblue"          "deepskyblue1"         "deepskyblue2"        
## [124] "deepskyblue3"         "deepskyblue4"         "dimgray"             
## [127] "dimgrey"              "dodgerblue"           "dodgerblue1"         
## [130] "dodgerblue2"          "dodgerblue3"          "dodgerblue4"         
## [133] "firebrick"            "firebrick1"           "firebrick2"          
## [136] "firebrick3"           "firebrick4"           "floralwhite"         
## [139] "forestgreen"          "gainsboro"            "ghostwhite"          
## [142] "gold"                 "gold1"                "gold2"               
## [145] "gold3"                "gold4"                "goldenrod"           
## [148] "goldenrod1"           "goldenrod2"           "goldenrod3"          
## [151] "goldenrod4"           "gray"                 "gray0"               
## [154] "gray1"                "gray2"                "gray3"               
## [157] "gray4"                "gray5"                "gray6"               
## [160] "gray7"                "gray8"                "gray9"               
## [163] "gray10"               "gray11"               "gray12"              
## [166] "gray13"               "gray14"               "gray15"              
## [169] "gray16"               "gray17"               "gray18"              
## [172] "gray19"               "gray20"               "gray21"              
## [175] "gray22"               "gray23"               "gray24"              
## [178] "gray25"               "gray26"               "gray27"              
## [181] "gray28"               "gray29"               "gray30"              
## [184] "gray31"               "gray32"               "gray33"              
## [187] "gray34"               "gray35"               "gray36"              
## [190] "gray37"               "gray38"               "gray39"              
## [193] "gray40"               "gray41"               "gray42"              
## [196] "gray43"               "gray44"               "gray45"              
## [199] "gray46"               "gray47"               "gray48"              
## [202] "gray49"               "gray50"               "gray51"              
## [205] "gray52"               "gray53"               "gray54"              
## [208] "gray55"               "gray56"               "gray57"              
## [211] "gray58"               "gray59"               "gray60"              
## [214] "gray61"               "gray62"               "gray63"              
## [217] "gray64"               "gray65"               "gray66"              
## [220] "gray67"               "gray68"               "gray69"              
## [223] "gray70"               "gray71"               "gray72"              
## [226] "gray73"               "gray74"               "gray75"              
## [229] "gray76"               "gray77"               "gray78"              
## [232] "gray79"               "gray80"               "gray81"              
## [235] "gray82"               "gray83"               "gray84"              
## [238] "gray85"               "gray86"               "gray87"              
## [241] "gray88"               "gray89"               "gray90"              
## [244] "gray91"               "gray92"               "gray93"              
## [247] "gray94"               "gray95"               "gray96"              
## [250] "gray97"               "gray98"               "gray99"              
## [253] "gray100"              "green"                "green1"              
## [256] "green2"               "green3"               "green4"              
## [259] "greenyellow"          "grey"                 "grey0"               
## [262] "grey1"                "grey2"                "grey3"               
## [265] "grey4"                "grey5"                "grey6"               
## [268] "grey7"                "grey8"                "grey9"               
## [271] "grey10"               "grey11"               "grey12"              
## [274] "grey13"               "grey14"               "grey15"              
## [277] "grey16"               "grey17"               "grey18"              
## [280] "grey19"               "grey20"               "grey21"              
## [283] "grey22"               "grey23"               "grey24"              
## [286] "grey25"               "grey26"               "grey27"              
## [289] "grey28"               "grey29"               "grey30"              
## [292] "grey31"               "grey32"               "grey33"              
## [295] "grey34"               "grey35"               "grey36"              
## [298] "grey37"               "grey38"               "grey39"              
## [301] "grey40"               "grey41"               "grey42"              
## [304] "grey43"               "grey44"               "grey45"              
## [307] "grey46"               "grey47"               "grey48"              
## [310] "grey49"               "grey50"               "grey51"              
## [313] "grey52"               "grey53"               "grey54"              
## [316] "grey55"               "grey56"               "grey57"              
## [319] "grey58"               "grey59"               "grey60"              
## [322] "grey61"               "grey62"               "grey63"              
## [325] "grey64"               "grey65"               "grey66"              
## [328] "grey67"               "grey68"               "grey69"              
## [331] "grey70"               "grey71"               "grey72"              
## [334] "grey73"               "grey74"               "grey75"              
## [337] "grey76"               "grey77"               "grey78"              
## [340] "grey79"               "grey80"               "grey81"              
## [343] "grey82"               "grey83"               "grey84"              
## [346] "grey85"               "grey86"               "grey87"              
## [349] "grey88"               "grey89"               "grey90"              
## [352] "grey91"               "grey92"               "grey93"              
## [355] "grey94"               "grey95"               "grey96"              
## [358] "grey97"               "grey98"               "grey99"              
## [361] "grey100"              "honeydew"             "honeydew1"           
## [364] "honeydew2"            "honeydew3"            "honeydew4"           
## [367] "hotpink"              "hotpink1"             "hotpink2"            
## [370] "hotpink3"             "hotpink4"             "indianred"           
## [373] "indianred1"           "indianred2"           "indianred3"          
## [376] "indianred4"           "ivory"                "ivory1"              
## [379] "ivory2"               "ivory3"               "ivory4"              
## [382] "khaki"                "khaki1"               "khaki2"              
## [385] "khaki3"               "khaki4"               "lavender"            
## [388] "lavenderblush"        "lavenderblush1"       "lavenderblush2"      
## [391] "lavenderblush3"       "lavenderblush4"       "lawngreen"           
## [394] "lemonchiffon"         "lemonchiffon1"        "lemonchiffon2"       
## [397] "lemonchiffon3"        "lemonchiffon4"        "lightblue"           
## [400] "lightblue1"           "lightblue2"           "lightblue3"          
## [403] "lightblue4"           "lightcoral"           "lightcyan"           
## [406] "lightcyan1"           "lightcyan2"           "lightcyan3"          
## [409] "lightcyan4"           "lightgoldenrod"       "lightgoldenrod1"     
## [412] "lightgoldenrod2"      "lightgoldenrod3"      "lightgoldenrod4"     
## [415] "lightgoldenrodyellow" "lightgray"            "lightgreen"          
## [418] "lightgrey"            "lightpink"            "lightpink1"          
## [421] "lightpink2"           "lightpink3"           "lightpink4"          
## [424] "lightsalmon"          "lightsalmon1"         "lightsalmon2"        
## [427] "lightsalmon3"         "lightsalmon4"         "lightseagreen"       
## [430] "lightskyblue"         "lightskyblue1"        "lightskyblue2"       
## [433] "lightskyblue3"        "lightskyblue4"        "lightslateblue"      
## [436] "lightslategray"       "lightslategrey"       "lightsteelblue"      
## [439] "lightsteelblue1"      "lightsteelblue2"      "lightsteelblue3"     
## [442] "lightsteelblue4"      "lightyellow"          "lightyellow1"        
## [445] "lightyellow2"         "lightyellow3"         "lightyellow4"        
## [448] "limegreen"            "linen"                "magenta"             
## [451] "magenta1"             "magenta2"             "magenta3"            
## [454] "magenta4"             "maroon"               "maroon1"             
## [457] "maroon2"              "maroon3"              "maroon4"             
## [460] "mediumaquamarine"     "mediumblue"           "mediumorchid"        
## [463] "mediumorchid1"        "mediumorchid2"        "mediumorchid3"       
## [466] "mediumorchid4"        "mediumpurple"         "mediumpurple1"       
## [469] "mediumpurple2"        "mediumpurple3"        "mediumpurple4"       
## [472] "mediumseagreen"       "mediumslateblue"      "mediumspringgreen"   
## [475] "mediumturquoise"      "mediumvioletred"      "midnightblue"        
## [478] "mintcream"            "mistyrose"            "mistyrose1"          
## [481] "mistyrose2"           "mistyrose3"           "mistyrose4"          
## [484] "moccasin"             "navajowhite"          "navajowhite1"        
## [487] "navajowhite2"         "navajowhite3"         "navajowhite4"        
## [490] "navy"                 "navyblue"             "oldlace"             
## [493] "olivedrab"            "olivedrab1"           "olivedrab2"          
## [496] "olivedrab3"           "olivedrab4"           "orange"              
## [499] "orange1"              "orange2"              "orange3"             
## [502] "orange4"              "orangered"            "orangered1"          
## [505] "orangered2"           "orangered3"           "orangered4"          
## [508] "orchid"               "orchid1"              "orchid2"             
## [511] "orchid3"              "orchid4"              "palegoldenrod"       
## [514] "palegreen"            "palegreen1"           "palegreen2"          
## [517] "palegreen3"           "palegreen4"           "paleturquoise"       
## [520] "paleturquoise1"       "paleturquoise2"       "paleturquoise3"      
## [523] "paleturquoise4"       "palevioletred"        "palevioletred1"      
## [526] "palevioletred2"       "palevioletred3"       "palevioletred4"      
## [529] "papayawhip"           "peachpuff"            "peachpuff1"          
## [532] "peachpuff2"           "peachpuff3"           "peachpuff4"          
## [535] "peru"                 "pink"                 "pink1"               
## [538] "pink2"                "pink3"                "pink4"               
## [541] "plum"                 "plum1"                "plum2"               
## [544] "plum3"                "plum4"                "powderblue"          
## [547] "purple"               "purple1"              "purple2"             
## [550] "purple3"              "purple4"              "red"                 
## [553] "red1"                 "red2"                 "red3"                
## [556] "red4"                 "rosybrown"            "rosybrown1"          
## [559] "rosybrown2"           "rosybrown3"           "rosybrown4"          
## [562] "royalblue"            "royalblue1"           "royalblue2"          
## [565] "royalblue3"           "royalblue4"           "saddlebrown"         
## [568] "salmon"               "salmon1"              "salmon2"             
## [571] "salmon3"              "salmon4"              "sandybrown"          
## [574] "seagreen"             "seagreen1"            "seagreen2"           
## [577] "seagreen3"            "seagreen4"            "seashell"            
## [580] "seashell1"            "seashell2"            "seashell3"           
## [583] "seashell4"            "sienna"               "sienna1"             
## [586] "sienna2"              "sienna3"              "sienna4"             
## [589] "skyblue"              "skyblue1"             "skyblue2"            
## [592] "skyblue3"             "skyblue4"             "slateblue"           
## [595] "slateblue1"           "slateblue2"           "slateblue3"          
## [598] "slateblue4"           "slategray"            "slategray1"          
## [601] "slategray2"           "slategray3"           "slategray4"          
## [604] "slategrey"            "snow"                 "snow1"               
## [607] "snow2"                "snow3"                "snow4"               
## [610] "springgreen"          "springgreen1"         "springgreen2"        
## [613] "springgreen3"         "springgreen4"         "steelblue"           
## [616] "steelblue1"           "steelblue2"           "steelblue3"          
## [619] "steelblue4"           "tan"                  "tan1"                
## [622] "tan2"                 "tan3"                 "tan4"                
## [625] "thistle"              "thistle1"             "thistle2"            
## [628] "thistle3"             "thistle4"             "tomato"              
## [631] "tomato1"              "tomato2"              "tomato3"             
## [634] "tomato4"              "turquoise"            "turquoise1"          
## [637] "turquoise2"           "turquoise3"           "turquoise4"          
## [640] "violet"               "violetred"            "violetred1"          
## [643] "violetred2"           "violetred3"           "violetred4"          
## [646] "wheat"                "wheat1"               "wheat2"              
## [649] "wheat3"               "wheat4"               "whitesmoke"          
## [652] "yellow"               "yellow1"              "yellow2"             
## [655] "yellow3"              "yellow4"              "yellowgreen"
colors()[179]  # View specific element of the color vector
## [1] "gray26"
colors()[179:190]  # View a selection of elements from the color vector
##  [1] "gray26" "gray27" "gray28" "gray29" "gray30" "gray31" "gray32"
##  [8] "gray33" "gray34" "gray35" "gray36" "gray37"

Another option: R Color Infographic

plot(x = data$year, y = data$realgdpgr, type = "p", col = colors()[145])  # or col='gold3'
plot(x = data$year, y = data$realgdpgr, type = "p", col = colors()[624])  # or col='tan4'

A Good Reference

plot(x = data$year, y = data$realgdpgr, type = "p", pch = 3)  # Change point style to crosses
plot(x = data$year, y = data$realgdpgr, type = "p", pch = 15)  # Change point style to filled squares
# Change point style to filled squares and increase point size to 3
plot(x = data$year, y = data$realgdpgr, type = "p", pch = 15, cex = 3)
plot(x = data$year, y = data$realgdpgr, type = "p", pch = "w")  # Change point style to 'w'
# Change point style to '$' and increase point size to 2
plot(x = data$year, y = data$realgdpgr, type = "p", pch = "$", cex = 2)
# Create a data.frame containing yearly average Real GDP Growth over all
# countries
library(plyr)
# Split-apply-combine via plyr
mean.growth <- ddply(.data = data, .variables = .(year), summarize, mean = mean(realgdpgr))
head(mean.growth)
##   year  mean
## 1 1960    NA
## 2 1961 5.612
## 3 1962 4.830
## 4 1963 5.310
## 5 1964 6.618
## 6 1965 5.010
# Line plot with solid line
plot(x = mean.growth$year, y = mean.growth$mean, type = "l", lty = 1)
# Line plot with medium dashed line
plot(x = mean.growth$year, y = mean.growth$mean, type = "l", lty = 2)
# Line plot with short dashed line
plot(x = mean.growth$year, y = mean.growth$mean, type = "l", lty = 3)
# Change line width to 2
plot(x = mean.growth$year, y = mean.growth$mean, type = "l", lty = 3, lwd = 2)
# Change line width to 3
plot(x = mean.growth$year, y = mean.growth$mean, type = "l", lty = 3, lwd = 3)
# Change line width to 4
plot(x = mean.growth$year, y = mean.growth$mean, type = "l", lty = 3, lwd = 4)

base/layering

# Subset data to create a few data.frames containing data for a single
# country
france.growth <- data[data$country == "France", ]
italy.growth <- data[data$country == "Italy", ]
spain.growth <- data[data$country == "Spain", ]
# First call to plot
plot(x = france.growth$year, y = france.growth$realgdpgr, type = "l", col = "red", 
    lwd = 2)
# Second call to plot
plot(x = italy.growth$year, y = italy.growth$realgdpgr, type = "l", col = "blue", 
    lwd = 2)
# First call to plot
plot(x = france.growth$year, y = france.growth$realgdpgr, type = "l", col = "red", 
    lwd = 2)
# First call to lines
lines(x = italy.growth$year, y = italy.growth$realgdpgr, type = "l", col = "blue", 
    lwd = 2)
# Second call to lines
lines(x = spain.growth$year, y = spain.growth$realgdpgr, type = "l", col = "darkgreen", 
    lwd = 2)
# First call to plot
plot(x = france.growth$year, y = france.growth$realgdpgr, type = "l", col = "red", 
    lwd = 2)
# First call to lines
lines(x = italy.growth$year, y = italy.growth$realgdpgr, type = "l", col = "blue", 
    lwd = 2)
# Second call to lines
lines(x = spain.growth$year, y = spain.growth$realgdpgr, type = "l", col = "darkgreen", 
    lwd = 2)
# First call to points
points(x = 1986, y = 6, pch = 13, col = colors()[116])
# Second call to points
points(x = 1986, y = -1, pch = 13, col = colors()[24])
# Third call to points
points(x = 1966, y = 2, pch = 13, col = colors()[391])
# Fourth call to points
points(x = 2008, y = 4, pch = 13, col = colors()[8])

base/options [annotations/reference lines/legends]

plot(x = france.growth$year, y = france.growth$realgdpgr, type = "l", col = "red", 
    lwd = 2)
lines(x = italy.growth$year, y = italy.growth$realgdpgr, type = "l", col = "blue", 
    lwd = 2)
lines(x = spain.growth$year, y = spain.growth$realgdpgr, type = "l", col = "darkgreen", 
    lwd = 2)
points(x = 1986, y = 6, pch = 13, col = colors()[116])
points(x = 1986, y = -1, pch = 13, col = colors()[24])
points(x = 1966, y = 2, pch = 13, col = colors()[391])
points(x = 2008, y = 4, pch = 13, col = colors()[8])
# First call to text
text(x = 1967, y = -1, labels = "No dot here")
plot(x = france.growth$year, y = france.growth$realgdpgr, type = "l", col = "red", 
    lwd = 2)
lines(x = italy.growth$year, y = italy.growth$realgdpgr, type = "l", col = "blue", 
    lwd = 2)
lines(x = spain.growth$year, y = spain.growth$realgdpgr, type = "l", col = "darkgreen", 
    lwd = 2)
points(x = 1986, y = 6, pch = 13, col = colors()[116])
points(x = 1986, y = -1, pch = 13, col = colors()[24])
points(x = 1966, y = 2, pch = 13, col = colors()[391])
points(x = 2008, y = 4, pch = 13, col = colors()[8])
text(x = 1967, y = -1, labels = "No dot here")
# First call to abline
abline(v = 1980, h = 0)
plot(x = france.growth$year, y = france.growth$realgdpgr, type = "l", col = "red", 
    lwd = 2)
lines(x = italy.growth$year, y = italy.growth$realgdpgr, type = "l", col = "blue", 
    lwd = 2)
lines(x = spain.growth$year, y = spain.growth$realgdpgr, type = "l", col = "darkgreen", 
    lwd = 2)
points(x = 1986, y = 6, pch = 13, col = colors()[116])
points(x = 1986, y = -1, pch = 13, col = colors()[24])
points(x = 1966, y = 2, pch = 13, col = colors()[391])
points(x = 2008, y = 4, pch = 13, col = colors()[8])
text(x = 1967, y = -1, labels = "No dot here")
abline(v = 1980, h = 0)
# First call to legend (note the vector position correspondence between each
# of the argument values)
legend("topright", legend = c("France", "Italy", "Spain"), col = c("red", "blue", 
    "darkgreen"), lty = c(1, 1, 1), lwd = c(2, 2, 2), cex = 0.8)

base/tables

par(mrow = c(ncols, nrows))
# STEP 1: Call 'par() for a 2x2 table'
par(mfrow = c(2, 2))

# STEP 2: Plot #1
plot(x = france.growth$year, y = france.growth$realgdpgr, type = "l", col = "red", 
    lwd = 2)
lines(x = italy.growth$year, y = italy.growth$realgdpgr, type = "l", col = "blue", 
    lwd = 2)
lines(x = spain.growth$year, y = spain.growth$realgdpgr, type = "l", col = "darkgreen", 
    lwd = 2)
points(x = 1986, y = 6, pch = 13, col = colors()[116])
points(x = 1986, y = -1, pch = 13, col = colors()[24])
points(x = 1966, y = 2, pch = 13, col = colors()[391])
points(x = 2008, y = 4, pch = 13, col = colors()[8])

# STEP 3: Plot #2
plot(x = france.growth$year, y = france.growth$realgdpgr, type = "l", col = "red", 
    lwd = 2)
lines(x = italy.growth$year, y = italy.growth$realgdpgr, type = "l", col = "blue", 
    lwd = 2)
lines(x = spain.growth$year, y = spain.growth$realgdpgr, type = "l", col = "darkgreen", 
    lwd = 2)
points(x = 1986, y = 6, pch = 13, col = colors()[116])
points(x = 1986, y = -1, pch = 13, col = colors()[24])
points(x = 1966, y = 2, pch = 13, col = colors()[391])
points(x = 2008, y = 4, pch = 13, col = colors()[8])
text(x = 1967, y = -1, labels = "No dot here")

# STEP 4: Plot #3
plot(x = france.growth$year, y = france.growth$realgdpgr, type = "l", col = "red", 
    lwd = 2)
lines(x = italy.growth$year, y = italy.growth$realgdpgr, type = "l", col = "blue", 
    lwd = 2)
lines(x = spain.growth$year, y = spain.growth$realgdpgr, type = "l", col = "darkgreen", 
    lwd = 2)
points(x = 1986, y = 6, pch = 13, col = colors()[116])
points(x = 1986, y = -1, pch = 13, col = colors()[24])
points(x = 1966, y = 2, pch = 13, col = colors()[391])
points(x = 2008, y = 4, pch = 13, col = colors()[8])
text(x = 1967, y = -1, labels = "No dot here")
abline(v = 1980, h = 0)

# STEP 5: Plot #4
plot(x = france.growth$year, y = france.growth$realgdpgr, type = "l", col = "red", 
    lwd = 2)
lines(x = italy.growth$year, y = italy.growth$realgdpgr, type = "l", col = "blue", 
    lwd = 2)
lines(x = spain.growth$year, y = spain.growth$realgdpgr, type = "l", col = "darkgreen", 
    lwd = 2)
points(x = 1986, y = 6, pch = 13, col = colors()[116])
points(x = 1986, y = -1, pch = 13, col = colors()[24])
points(x = 1966, y = 2, pch = 13, col = colors()[391])
points(x = 2008, y = 4, pch = 13, col = colors()[8])
text(x = 1967, y = -1, labels = "No dot here")
abline(v = 1980, h = 0)
legend("topright", legend = c("France", "Italy", "Spain"), col = c("red", "blue", 
    "darkgreen"), lty = c(1, 1, 1), lwd = c(2, 2, 2), cex = 0.8)

Out with the old...

lattice v. ggplot2

  1. faster (though only noticeable over many and large plots)

  2. simpler (at first)

  3. better at trellis graphs

  4. able to do 3d graphs

  1. generally more elegant

  2. more syntactically logical (and therefore simpler, once you learn it)

  3. better at grouping

  4. able to interface with maps

Basic usage: lattice

The general call for lattice graphics looks something like this:

graph_type(formula, data=, [options])

The specifics of the formula differ for each graph type, but the general format is straightforward

y  # Show the distribution of y

y ~ x  # Show the relationship between x and y 

y ~ x | A  # Show the relationship between x and y conditional on the values of A

y ~ x | A * B  # Show the relationship between x and y conditional on the combinations of A and B

z ~ y * x  # Show the 3D relationship between x, y, and z

Basic usage: ggplot2

The general call for ggplot2 graphics looks something like this:

ggplot(data=, aes(x=,y=, [options]))+geom_xxxx()+...+...+...

Note that ggplot2 graphs in layers in a continuing call (hence the endless +...+...+...), which really makes the extra layer part of the call

...+geom_xxxx(data=, aes(x=,y=,[options]),[options])+...+...+...

You can see the layering effect by comparing the same graph with different colors for each layer

ggplot(data = data, aes(x = vturn, y = realgdpgr)) + geom_point(color = "black") + 
    geom_point(aes(x = vturn, y = unemp), color = "red")
ggplot(data = data, aes(x = vturn, y = realgdpgr)) + geom_point(color = "red") + 
    geom_point(aes(x = vturn, y = unemp), color = "black")

Comparing lattice and ggplot

lattice v. ggplot2: Densities

densityplot(~vturn, data = data)  # lattice
ggplot(data = data, aes(x = vturn)) + geom_density()  # ggplot2

lattice v. ggplot2: X-Y scatter plots

xyplot(outlays ~ vturn, data = data)  # lattice
ggplot(data = data, aes(x = vturn, y = outlays)) + geom_point()  # ggplot2

lattice v. ggplot2: X-Y line plots

xyplot(mean ~ year, data = mean.growth, type = "l")  # lattice
ggplot(data = mean.growth, aes(x = year, y = mean)) + geom_line()  # ggplot2 

lattice v. ggplot2: bar plots

# Create data.frame of average growth rates by country over time
growth <- ddply(.data = data, .variables = .(country), summarize, mean = mean(realgdpgr, 
    na.rm = T))

barchart(mean ~ country, data = growth)  # lattice
ggplot(data = growth, aes(x = country, y = mean)) + geom_bar()  # ggplot2

lattice v. ggplot2: box plots

bwplot(vturn ~ country, data = data)  # lattice
ggplot(data = data, aes(x = country, y = vturn)) + geom_boxplot()  # ggplot2

lattice v. ggplot2: "trellis" plots

xyplot(vturn ~ year | country, data = data)  # lattice
ggplot(data = data, aes(x = year, y = vturn)) + geom_point() + facet_wrap(~country)  # ggplot2

lattice v. ggplot2: countour plots

data(volcano)  # Load volcano contour data
volcano[1:10, 1:10]  # Examine volcano dataset (first 10 rows and columns)
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
##  [1,]  100  100  101  101  101  101  101  100  100   100
##  [2,]  101  101  102  102  102  102  102  101  101   101
##  [3,]  102  102  103  103  103  103  103  102  102   102
##  [4,]  103  103  104  104  104  104  104  103  103   103
##  [5,]  104  104  105  105  105  105  105  104  104   103
##  [6,]  105  105  105  106  106  106  106  105  105   104
##  [7,]  105  106  106  107  107  107  107  106  106   105
##  [8,]  106  107  107  108  108  108  108  107  107   106
##  [9,]  107  108  108  109  109  109  109  108  108   107
## [10,]  108  109  109  110  110  110  110  109  109   108
volcano3d <- melt(volcano)  # Use reshape2 package to melt the data
head(volcano3d)  # Examine volcano3d dataset (head)
##   Var1 Var2 value
## 1    1    1   100
## 2    2    1   101
## 3    3    1   102
## 4    4    1   103
## 5    5    1   104
## 6    6    1   105
names(volcano3d) <- c("xvar", "yvar", "zvar")  # Rename volcano3d columns
contourplot(zvar ~ xvar + yvar, data = volcano3d)  # lattice
ggplot(data = volcano3d, aes(x = xvar, y = yvar, z = zvar)) + geom_contour()  # ggplot2

lattice v. ggplot2: tile/image/level plots

levelplot(zvar ~ xvar + yvar, data = volcano3d)  # lattice
ggplot(data = volcano3d, aes(x = xvar, y = yvar, z = zvar)) + geom_tile(aes(fill = zvar))  # ggplot2

lattice: 3D plots

# Create a subset of the dataset containing only data for France
france.data <- data[data$country == "France", ]
cloud(outlays ~ year * vturn, data = france.data)
# Create a subset of the dataset containing only data for Greece, Portugal,
# Ireland, and Spain
pigs.data <- data[data$country %in% c("Greece", "Portugal", "Ireland", "Spain"), 
    ]
cloud(outlays ~ year * vturn | country, data = pigs.data)

ggplot2: Panel plots

ggplot(data = pigs.data, aes(x = year, y = vturn, color = country)) + geom_line()

lattice v. ggplot2: options [labeling]

xyplot(outlays ~ vturn, data = data, xlab = "Voter Turnout (%)", ylab = "Government Outlays", 
    main = "This Graph is Also Great")  # lattice
ggplot(data = data, aes(x = vturn, y = outlays)) + geom_point() + xlab(label = "Voter Turnout (%)") + 
    ylab(label = "Government Outlays") + ggtitle(label = "This Graph is Also Also Great")  # ggplot2

lattice v. ggplot2: options [axis + size scaling]

xyplot(outlays ~ vturn, data = data, xlim = c(80, 100))  # lattice
xyplot(outlays ~ vturn, data = data, xlim = c(80, 100), cex = 2)  # lattice
xyplot(outlays ~ vturn, data = data, xlim = c(80, 100), cex = 0.5)  # lattice
ggplot(data = data, aes(x = vturn, y = outlays)) + geom_point() + xlim(80, 100)  # ggplot2
ggplot(data = data, aes(x = vturn, y = outlays)) + geom_point(size = 3) + xlim(80, 
    100)  # ggplot2
ggplot(data = data, aes(x = vturn, y = outlays)) + geom_point(size = 1) + xlim(80, 
    100)  # ggplot2

lattice v. ggplot2: options [graphical parameters]

xyplot(outlays ~ vturn, data = data, col = colors()[145])  #lattice
xyplot(outlays ~ vturn, data = data, col = "red")  #lattice
ggplot(data = data, aes(x = vturn, y = outlays)) + geom_point(color = colors()[145])  # ggplot2
ggplot(data = data, aes(x = vturn, y = outlays)) + geom_point(color = "red")  # ggplot2
xyplot(outlays ~ vturn, data = data, pch = 3)  # lattice
xyplot(outlays ~ vturn, data = data, pch = 15)  # lattice
ggplot(data = data, aes(x = vturn, y = outlays)) + geom_point(shape = 3)  # ggplot2
ggplot(data = data, aes(x = vturn, y = outlays)) + geom_point(shape = 15)  # ggplot2
xyplot(outlays ~ vturn, data = data, pch = 3)  # lattice
xyplot(outlays ~ vturn, data = data, pch = 15)  # lattice
xyplot(outlays ~ vturn, data = data, pch = "w")  # lattice
xyplot(outlays ~ vturn, data = data, pch = "$", cex = 2)  # lattice
ggplot(data = data, aes(x = vturn, y = outlays)) + geom_point(shape = 3)  # ggplot2
ggplot(data = data, aes(x = vturn, y = outlays)) + geom_point(shape = 15)  # ggplot2
ggplot(data = data, aes(x = vturn, y = outlays)) + geom_point(shape = "w")  # ggplot2
ggplot(data = data, aes(x = vturn, y = outlays)) + geom_point(shape = "$", size = 5)  # ggplot2
xyplot(outlays ~ vturn, data = data, type = "l", lty = 1)  # lattice
xyplot(outlays ~ vturn, data = data, type = "l", lty = 2)  # lattice
xyplot(outlays ~ vturn, data = data, type = "l", lty = 3)  # lattice
xyplot(outlays ~ vturn, data = data, type = "l", lty = 3, lwd = 2)  # lattice
xyplot(outlays ~ vturn, data = data, type = "l", lty = 3, lwd = 3)  # lattice
xyplot(outlays ~ vturn, data = data, type = "l", lty = 3, lwd = 4)  # lattice
ggplot(data = data, aes(x = vturn, y = outlays)) + geom_line(linetype = 1)  # ggplot2
ggplot(data = data, aes(x = vturn, y = outlays)) + geom_line(linetype = 2)  # ggplot2
ggplot(data = data, aes(x = vturn, y = outlays)) + geom_line(linetype = 3)  # ggplot2
ggplot(data = data, aes(x = vturn, y = outlays)) + geom_line(linetype = 3, size = 1)  # ggplot2
ggplot(data = data, aes(x = vturn, y = outlays)) + geom_line(linetype = 3, size = 1.5)  # ggplot2
ggplot(data = data, aes(x = vturn, y = outlays)) + geom_line(linetype = 3, size = 2)  # ggplot2

ggplot2 and the Grammar of Graphics

ggplot2 and the Grammar of Graphics

  1. One or more statistics conveying information about the data (identities, means, medians, etc.)

  2. A coordinate system that differentiates between the intersections of statistics (at most two for ggplot, three for lattice)

  3. Geometries that differentiate between off-coordinate variation in kind

  4. Scales that differentiate between off-coordinate variation in degree

Anatomy of aes()

ggplot(data = , aes(x = , y = , color = , linetype = , shape = , size = ))

ggplot2 is optimized for showing variation on all four aesthetic types

# Differences in kind using color
ggplot(data = pigs.data, aes(x = year, y = vturn)) + geom_line(aes(color = country))

Note what happens when we specify the color parameter outside of the aesthetic operator. ggplot2 views these specifications as invalid graphical parameters.

ggplot(data = pigs.data, aes(x = year, y = vturn)) + geom_line(color = country)
## Error: object 'country' not found
ggplot(data = pigs.data, aes(x = year, y = vturn)) + geom_line(color = "country")
## Error: invalid color name 'country'
ggplot(data = pigs.data, aes(x = year, y = vturn)) + geom_line(color = "red")
# Differences in kind using line types
ggplot(data = pigs.data, aes(x = year, y = vturn)) + geom_line(aes(linetype = country))
# Differences in kind using point shapes
ggplot(data = pigs.data, aes(x = year, y = vturn)) + geom_point(aes(shape = country))
# Differences in degree using color
ggplot(data = pigs.data, aes(x = year, y = realgdpgr)) + geom_point(aes(color = vturn))
# Differences in degree using point size
ggplot(data = pigs.data, aes(x = year, y = realgdpgr)) + geom_point(aes(size = vturn))
# Multiple non-cartesian aesthetics (differences in kind using color, degree
# using point size)
ggplot(data = pigs.data, aes(x = year, y = realgdpgr)) + geom_point(aes(color = country, 
    size = vturn))

Fitted lines and curves with ggplot2

ggplot(data = pigs.data, aes(x = year, y = vturn)) + geom_point()
# Add linear model (lm) smoother
ggplot(data = pigs.data, aes(x = year, y = vturn)) + geom_point() + geom_smooth(method = "lm")
# Add local linear model (loess) smoother, span of 0.75
ggplot(data = pigs.data, aes(x = year, y = vturn)) + geom_point() + geom_smooth(method = "loess", 
    span = 0.75)
# Add local linear model (loess) smoother, span of 0.25
ggplot(data = pigs.data, aes(x = year, y = vturn)) + geom_point() + geom_smooth(method = "loess", 
    span = 0.25)
# Add linear model (lm) smoother, no standard error shading
ggplot(data = pigs.data, aes(x = year, y = vturn)) + geom_point() + geom_smooth(method = "lm", 
    se = F)
# Add local linear model (loess) smoother, no standard error shading
ggplot(data = pigs.data, aes(x = year, y = vturn)) + geom_point() + geom_smooth(method = "loess", 
    se = F)
# Add a local linear (loess) smoother for each country
ggplot(data = pigs.data, aes(x = year, y = vturn)) + geom_point(aes(color = country)) + 
    geom_smooth(aes(color = country))
# Add a local linear (loess) smoother for each country, no standard error
# shading
ggplot(data = pigs.data, aes(x = year, y = realgdpgr)) + geom_point(aes(color = country, 
    size = vturn)) + geom_smooth(aes(color = country), se = F)

lattice v. ggplot2: tables

# Initialize gridExtra library
library(gridExtra)
# Create 3 plots to combine in a table
plot1 <- ggplot(data = pigs.data, aes(x = year, y = vturn, color = )) + geom_line(aes(color = country))
plot2 <- ggplot(data = pigs.data, aes(x = year, y = vturn, linetype = )) + geom_line(aes(linetype = country))
plot3 <- ggplot(data = pigs.data, aes(x = year, y = vturn, shape = )) + geom_point(aes(shape = country))
# Call grid.arrange
grid.arrange(plot1, plot2, plot3, nrow = 3, ncol = 1)

Exporting

Two basic image types

  1. Raster/Bitmap (.png, .jpeg)

Every pixel of a plot contains its own separate coding; not so great if you want to resize the image

jpeg(filename = "example.png", width = , height = )
plot(x, y)
dev.off()
  1. Vector (.pdf, .ps)

Every element of a plot is encoded with a function that gives its coding conditional on several factors; great for resizing

pdf(filename = "example.pdf", width = , height = )
plot(x, y)
dev.off()

Exporting with lattice v. ggplot

# Assume we saved our plot is an object called example.plot

# lattice
trellis.device(device = "pdf", filename = "example.pdf")
print(example.plot)
dev.off()

# ggplot2
ggsave(filename = "example.pdf", plot = example.plot, scale = , width = , height = )  # ggplot2

stop()